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JEFFREY R. CHASNOV 


DIFFERENTIAL 
EQUATIONS 


FOR ENGINEERS 


Preface 


What follows are my lecture notes for a first course in differential equations, taught 
at the Hong Kong University of Science and Technology. Included in these notes 
are links to short tutorial videos posted on YouTube. 

Much of the material of Chapters 2-6 and 8 has been adapted from the widely 
used textbook “Elementary differential equations and boundary value problems” 
by Boyce & DiPrima (John Wiley & Sons, Inc., Seventh Edition, ©2001). Many of 
the examples presented in these notes may be found in this book. The material of 
Chapter 7 is adapted from the textbook “Nonlinear dynamics and chaos” by Steven 
H. Strogatz (Perseus Publishing, ©1994). 


All web surfers are welcome to download these notes, watch the YouTube videos, 
and to use the notes and videos freely for teaching and learning. 


I also have some online courses on Coursera. A lot of time and effort has gone 
into their production, and the video lectures have better video quality than the ones 
prepared for these notes. You can click on the links below to explore these courses. 
If you want to learn differential equations, have a look at 

Differential Equations for Engineers 
If your interests are matrices and elementary linear algebra, try 


Matrix Algebra for Engineers 


If you want to learn vector calculus (also known as multivariable calculus, or calcu- 
lus three), you can sign up for 


Vector Calculus for Engineers 
And if your interest is numerical methods, have a go at 


Numerical Methods for Engineers 


JEFFREY R. CHASNOV 
Hong Kong 
February 2021 
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Chapter 0 


A short mathematical review 


A basic understanding of calculus is required to undertake a study of differential 
equations. This zero chapter presents a short review. 


0.1 The trigonometric functions 
The Pythagorean trigonometric identity is 

sin? x + cos* x = 1, 
and the addition theorems are 


sin(x + y) = sin(x) cos(y) + cos(x) sin(y), 
cos(x + y) = cos(x) cos(y) — sin(x) sin(y). 


Also, the values of sinx in the first quadrant can be remembered by the rule of 
quarters, with 0° = 0, 30° = 77/6, 45° = 71/4, 60° = 77/3, 90° = 7/2: 


a test. OE noes 342 
sin 0 -/% sin 30° = 1 sin 45° = 1 


4 
sin 60° = . sin 90° = re 
The following symmetry properties are also useful: 
sin(77/2—x) =cosx, cos(7/2— x) =sinx; 


and 
sin(—x) = —sin(x), cos(—x) = cos(x). 


0.2 The exponential function and the natural logarithm 


The transcendental number e, approximately 2.71828, is defined as 


1 n 
e= lim (1 + x) : 
n—- oo n 


The exponential function exp (x) = e* and natural logarithm In x are inverse func- 
tions satisfying 

e"*=x, Ine* =x. 
The usual rules of exponents apply: 


ee = ety, e*/e! — ey, (e*)P =— eP*. 
The corresponding rules for the logarithmic function are 


In(xy) =Inx+Iny, In(x/y)=Inx—Iny, Inx? = pInx. 


0.3. DEFINITION OF THE DERIVATIVE 


0.3 Definition of the derivative 


The derivative of the function y = f(x), denoted as f(x) or dy/dx, is defined as 
the slope of the tangent line to the curve y = f(x) at the point (x,y). This slope is 
obtained by a limit, and is defined as 


h-0 h 


0.4 Differentiating a combination of functions 


0.4.1 The sum or difference rule 
The derivative of the sum of f(x) and g(x) is 
(f+s) =f +8’. 


Similarly, the derivative of the difference is 
Ug =f 8: 


0.4.2 The product rule 
The derivative of the product of f(x) and g(x) is 
(fg) = fist fe’, 
and should be memorized as “the derivative of the first times the second plus the 
first times the derivative of the second.” 
0.4.3 The quotient rule 
The derivative of the quotient of f(x) and g(x) is 


(fy = gt 


y 


and should be memorized as “the derivative of the top times the bottom minus the 
top times the derivative of the bottom over the bottom squared.” 


0.4.4 The chain rule 


The derivative of the composition of f(x) and g(x) is 


and should be memorized as “the derivative of the outside times the derivative of 
the inside.” 
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0.5. DIFFERENTIATING ELEMENTARY FUNCTIONS 


0.5 Differentiating elementary functions 


0.5.1 The power rule 


The derivative of a power of x is given by 


d 
an = px, 


0.5.2 Trigonometric functions 
The derivatives of sin x and cos x are 
(sinx)' =cosx, (cosx)’ = —sinx. 


We thus say that “the derivative of sine is cosine,” and “the derivative of cosine is 
minus sine.” Notice that the second derivatives satisfy 


(sinx)” = —sinx, (cosx)” = —cosx. 


0.5.3 Exponential and natural logarithm functions 


The derivative of e* and Inx are 


1 
x\t x / sae 
(e*)' =e*, (Inx) = 


0.6 Definition of the integral 


The definite integral of a function f(x) > 0 from x = a to b (b > a) is defined 
as the area bounded by the vertical lines x = a, x = b, the x-axis and the curve 
y = f(x). This “area under the curve” is obtained by a limit. First, the area is 
approximated by a sum of rectangle areas. Second, the integral is defined to be the 
limit of the rectangle areas as the width of each individual rectangle goes to zero 
and the number of rectangles goes to infinity. This resulting infinite sum is called a 
Riemann Sum, and we define 


N 
[flax = tim Y flat (n—1)M) +h, (2) 
a 30 (4 


where N = (b —a)/h is the number of terms in the sum. The symbols on the left- 
hand-side of (2) are read as “the integral from a to b of f of x dee x.” The Riemann 
Sum definition is extended to all values of a and b and for all values of f (x) (positive 
and negative). Accordingly, 


[ feiax =- [soar and [rene = [ seoae 
Also, 
[ feax= | f(xdx + [ f()dx, 


which states when f(x) > 0 anda < b <c that the total area is equal to the sum of 
its parts. 
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0.7. THE FUNDAMENTAL THEOREM OF CALCULUS 


0.7. The fundamental theorem of calculus 
View tutorial on YouTube 


Using the definition of the derivative, we differentiate the following integral: 


2" f(s)ds — [* f(s)ds 
h 


i f(s)ds 
h 


This result is called the fundamental theorem of calculus, and provides a connection 
between differentiation and integration. 

The fundamental theorem teaches us how to integrate functions. Let F(x) be a 
function such that F’(x) = f(x). We say that F(x) is an antiderivative of f(x). Then 
from the fundamental theorem and the fact that the derivative of a constant equals 
Zero, 


F(x) = [ fos +e. 


Now, F(a) = c and F(b) = f[ : f(s)ds + F(a). Therefore, the fundamental theorem 
shows us how to integrate a function f(x) provided we can find its antiderivative: 


[ flsvas =F) - Fla). @) 


Unfortunately, finding antiderivatives is much harder than finding derivatives, and 
indeed, most complicated functions cannot be integrated analytically. 

We can also derive the very important result (3) directly from the definition of 
the derivative (1) and the definite integral (2). We will see it is convenient to choose 
the same hi in both limits. With F’(x) = f(x), we have 


n=1 
N. F(a+nh) — F(a+ (n—1)h) 
=f h i 
ON 
= Le + nh) — F(a + (n—1)h) 


The last expression has an interesting structure. All the values of F(x) evaluated 
at the points lying between the endpoints a and b cancel each other in consecutive 
terms. Only the value —F(a) survives when n = 1, and the value +F(b) when 
n = N, yielding again (3). 
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0.8. DEFINITE AND INDEFINITE INTEGRALS 


0.8 Definite and indefinite integrals 


The Riemann sum definition of an integral is called a definite integral. It is convenient 
to also define an indefinite integral by 


| flax = F(x), 


where F(x) is the antiderivative of f(x). 


0.9 Indefinite integrals of elementary functions 


From our known derivatives of elementary functions, we can determine some sim- 
ple indefinite integrals. The power rule gives us 


ntl 
n —_ = 
[x= +0, n#-—l. 


When n = —1, and x is positive, we have 


[cars tnx te. 
x 


If x is negative, using the chain rule we have 


d 1 
aes In ( — x) = 6 
Therefore, since 


i= —x ifx <0; 
~ |) x ifx>0, 


we can generalize our indefinite integral to strictly positive or strictly negative x: 
1 
/ —dx =In|x| +c. 
x 
Trigonometric functions can also be integrated: 
J cosxdx =sinx+c, |[sinxax =—cosx+c. 
Easily proved identities are an addition rule: 
[UF + s@o)ax= f fede [ goaz 
and multiplication by a constant: 
[ Af(aax = A | f(x)dx. 
This permits integration of functions such as 
x? 7x? 


[2+7x+2)ax= 24 +2x+c, 


[cos + sinx)dx = 5sinx —cosx+c. 
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0.10. SUBSTITUTION 


0.10 Substitution 


More complicated functions can be integrated using the chain rule. Since 


# f(g(x)) = f(g) -8'), 


we have 
ea (g(x)) - 9’ (x)dx = f(g(x)) +e. 


This integration formula is usually implemented by letting y = g(x). Then one 
writes dy = 9/(x)dx to obtain 
[fF (ed)! eoax= fF wdy 
= fly) +e 
= f(g(x)) +e. 


0.11 Integration by parts 


Another integration technique makes use of the product rule for differentiation. 
Since 


(fg) = figt fe, 
we have 

f'g = (fg) — fs’. 
Therefore, 

| fF @)s@)dx = feds) - | F)g!(Xdx. 
Commonly, the above integral is done by writing 
uw = g(x) dv = f'(x)dx 
du=g'(x)dx v=f(x). 


Then, the formula to be memorized is 


[uae =uv— [eau. 


0.12 ‘Taylor series 


A Taylor series of a function f(x) about a point x = a is a power series repre- 
sentation of f(x) developed so that all the derivatives of f(x) at a match all the 
derivatives of the power series. Without worrying about convergence here, we have 


f(x) = f(a) + f’(a)(x —a) 4 ie (x —a)?4 ie (x—a)o +... 


Notice that the first term in the power series matches f(a), all other terms vanishing, 
the second term matches /'(a), all other terms vanishing, etc. Commonly, the Taylor 
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0.13. FUNCTIONS OF SEVERAL VARIABLES 


series is developed with a = 0. We will also make use of the Taylor series in a 
slightly different form, with x = x, +¢€ and a= x;,: 


flor te) = flr) + fimer Le, Lay 


Another way to view this series is that of g(€) = f(x. +¢€), expanded about e = 0. 
Taylor series that are commonly used include 


2 3 
ia x x 
e SE op tap tke 
a xe x? 
SU Se ag ep Ne 
2 x2 x4 
cosx = ae ae 
1 
Teietie aS np. Tor | ye; 
2 3 
In(l+x)=x 5 tS ..., for |x| <1. 


0.13 Functions of several variables 


For simplicity, we consider a function f = f(x,y) of two variables, though the 
results are easily generalized. The partial derivative of f with respect to x is defined 


as 
OF tn SEH) — Flay) 
ox h-0 h 


and similarly for the partial derivative of f with respect to y. To take the partial 
derivative of f with respect to x, say, take the derivative of f with respect to x 
holding y fixed. As an example, consider 


f(xy) =2ey +y¥°. 


We have 


a = 6x"y’, - = 4x? y + By’. 


ox 


Second derivatives are defined as the derivatives of the first derivatives, so we have 
of 2 Of 3 
aya — lax: aed + 6Y; 


and the mixed second partial derivatives are 


of Re A 2 
Soy 12x*y, yee 12x*y. 


In general, mixed partial derivatives are independent of the order in which the 
derivatives are taken. 
Partial derivatives are necessary for applying the chain rule. Consider 


df = f(x +dx,y + dy) — f(x,y). 
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0.14. COMPLEX NUMBERS 


We can write df as 


df = [f(x +dx,y + dy) — f(x,y + dy)] + [f(xy + dy) — f(x,y)] 


J pe 
- 54% + ay 


If one has f = f(x(t),y(t)), say, then 


df _ofdx | of dy 
dt  axdt ' dy dt’ 


And if one has f = f(x(r,@),y(r,@)), say, then 


af afax , afay af _ afax , afay 
or oxor  dyor’ 00 axdd dyad 


A Taylor series of a function of several variables can also be developed. Here, all 
partial derivatives of f(x,y) at (a,b) match all the partial derivatives of the power 
series. With the notation 


_ of _ of Of _ &f _ oF 
fx 7 et fy _ oy’ fx = ae" fry = axay’ fyy — ay?’ etc., 


we have 


f(x,y) = f (a,b) + fr(a,b)(x — a) + fy (a,b) (y — b) 
+ x (fax (a,b) (x = a)* + 2fry(a,b) (x —a)(y —b) + fyy(4, b)(y— b)?) he 


0.14 Complex numbers 


View tutorial on YouTube: Complex Numbers 
View tutorial on YouTube: Complex Exponential Function 


We define the imaginary number i to be one of the two numbers that satisfies the 
rule (i ir = —1, the other number being —i. Formally, we write i = VJ—-1.A complex 
number z is written as 

z=x+y, 


where x and y are real numbers. We call x the real part of z and y the imaginary 
part and write 
x=Rez, y=Imz. 


Two complex numbers are equal if and only if their real and imaginary parts are 
equal. 
The complex conjugate of z = x + iy, denoted as Z, is defined as 
Z=x-—ly. 
Using z and Z, we have 


Rez = 5 (z+2), Imz = = (z—2). (4) 
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0.14. COMPLEX NUMBERS 


Furthermore, 


zz = (x +iy)(x — iy) 
eee Py? 
eee iP; 


and we define the absolute value of z, also called the modulus of z, by 


|z| = (z2)'/? 


= ,/xr24y2. 


We can add, subtract, multiply and divide complex numbers to get new complex 
numbers. With z = x + iy and w =s + it, and x,y,s,t real numbers, we have 


Furthermore, 


and 


Similarly 


s)+i(y+t); z-—w=(x-s)+i(y—t); 


zw = (x+iy)(s + it) 
= (xs — yt) + i(xt+ ys); 


Z zw 
w wo 
_ (x +iy)(s — it) 
gee 
(xe ryt) YS = 4E) 
= eae gape 


|zw| = Vs — yt)? + (xt + ys)? 
= +P)(e+P) 


= |z||wl; 


ZW = (xs — yt) —i(xt+ ys) 
(x — iy)(s — it) 


ZW. 


pom Gas 


Also, Z-- w = Z+W. However, |z + w| < |z| + |w|, a theorem known as the triangle 


inequality. 


It is especially interesting and useful to consider the exponential function of an 
imaginary argument. Using the Taylor series of an exponential function, we have 


e? — 1+ (i0) 4 


(19)? (i0)> , (i0)* | (i0)5 
2° Bho AEB 


ge Oe , een 
=( TRAE es -i(6 att 5 ba) 


= cos@+isiné. 
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0.14. COMPLEX NUMBERS 


Since we have determined that 
cos@ = Ree”, sin@ = Ime”, (5) 


we also have using (4) and (5), the frequently used expressions 


ei? e710 i? _ ,—i0 
cos 6 = — sin @ = oe 
The much celebrated Euler’s identity derives from e = cos +isind by setting 
6 = 7, and using cos 7 = —1 and sinz = 0: 
7 411=0, 


and this identity links the five fundamental numbers—0, 1, i, e and 7—using three 
basic mathematical operations—addition, multiplication and exponentiation—only 
once. 


y Z=X+ly 


Im(z) 


Figure 1: The complex plane. 


The complex number z can be represented in the complex plane with Rez as the 
x-axis and Imz as the y-axis (see Fig. 1). This leads to the polar representation of 
z=x+iy: 

z= rel”, 
where r = |z| and tan@ = y/x. We define argz = 0. Note that 6 is not unique, 
though it is conventional to choose the value such that —7 < 6 < mz, and @=0 
when r = 0. 

The polar form of a complex number can be useful when multiplying numbers. 
For example, if z; = ret and Z) = rete, then 21Z2 = ry reli t 02) | In particular, if 
’2 = 1, then multiplication of z; by z2 spins the representation of z; in the complex 
plane an angle @2 counterclockwise. 

Useful trigonometric relations can be derived using e and properties of the 
exponential function. The addition law can be derived from 


eilxty) = eX ely. 
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0.14. COMPLEX NUMBERS 


We have 
cos(x + y) t+isin(x + y) = (cosx +isinx)(cosy + isiny) 
= (cos x cosy — sinxsiny) + i(sinx cosy + cos x sin y); 
yielding 
cos(x + y) =cosxcosy—sinxsiny, sin(x+y) =sinxcosy+cosxsiny. 
De Moivre’s Theorem derives from e'”* = (e'®)", yielding the identity 
cos(n@) + isin(n@) = (cos@+isin@)". 
For example, if n = 2, we derive 
cos 26 + isin 20 = (cos@ +isin@)? 
= (cos* @ — sin? @) + 2icos @siné. 
Therefore, 


cos 20 = cos* 6 — sin? 6, sin20 = 2cos@sin0. 


Example: Write i as a standard complex number 


To solve this example, we first need to define what is meant by the square root 
of a complex number. The meaning of \/z is the complex number whose square 
is z. There will always be two such numbers, because (\/z)* = (—/z)? = z. One 
can not define the positive square root because complex numbers are not defined 
as positive or negative. 

We will show two methods to solve this problem. The first most straightforward 
method writes 


Vi=x+ iy. 
Squaring both sides, we obtain 
i= x? = + 2xyi; 


and equating the real and imaginary parts of this equation yields the two real equa- 
tions 
x7-y =0, 2xy =1. 


The first equation yields y = +x. With y = x, the second equation yields 2x2 = 1 
with two solutions x = £V2/2. With y = —x, the second equation yields 2x2 = 1, 
which has no solution for real x. We have therefore found that 


. V2 on 2 
vi=+ (+2). 


The second solution method makes use of the polar form of complex numbers. 
The algebra required for this method is somewhat simpler, especially for finding 
cube roots, fourth roots, etc. We know that i = elTt/ 2 but more generally because of 
the periodic nature of the polar angle, we can write 


i= ell Z+2mk) 
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where k is an integer. We then have 


Vi = ji/2 = el 447k) = eit pit/4 = +eit/4 


where we have made use of the usual properties of the exponential function, and 


eitk — +] for k even or odd. Converting back to standard form, we have 


Vi=t (cos 7/4+isin7/4) = + & +i) ; 


The fundamental theorem of algebra states that every polynomial equation of 
degree n has exactly n complex roots, counted with multiplicity. Two familiar ex- 
amples would be x* — 1 = (x + 1)(x — 1) = 0, with two roots x, = —1 and x2 = 1; 
and x? — 2x +1 = (x —1)* = 0, with one root x; = 1 with multiplicity two. 

The problem of finding the nth roots of unity is to solve the polynomial equation 


Zz? =1 


for the n complex values of z. We have z; = 1 for n = 1; and z; = 1, z2 = —1 for 
n = 2. Beyond n = 2, some of the roots are complex and here we find the cube 
roots of unity, that is, the three values of z that satisfy P=. Writing 1 = el27tk 
where k is an integer, we have 


1; 


1 
ee 27/3. 
y 


—_ ei2nk/3 = e! 


z=(1)3= (eo) 
elArt/3_ 
Using cos (27/3) = —1/2, sin (27/3) = V3/2, cos (47/3) = —1/2, sin (47¢/3) = 
—,/3/2, the three cube roots of unity are given by 


1 V3 1 V3 


z=1, ma —5 tis 23 = -= -—i—— 


These three roots are evenly spaced around the unit circle in the complex plane, as 
shown in the figure below. 


ei2n/3 


ei4n/3 
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Chapter 1 


Introduction to odes 


A differential equation is an equation for a function that relates the values of the 
function to the values of its derivatives. An ordinary differential equation (ode) is a 
differential equation for a function of a single variable, e.g., x(t), while a partial dif- 
ferential equation (pde) is a differential equation for a function of several variables, 
e.g., v(x, y,z,t). An ode contains ordinary derivatives and a pde contains partial 
derivatives. Typically, pde’s are much harder to solve than ode’s. 


1.1 The simplest type of differential equation 
View tutorial on YouTube 


The simplest ordinary differential equations can be integrated directly by finding 
antiderivatives. These simplest odes have the form 


d"x 
dt" Ge) 


where the derivative of x = x(t) can be of any order, and the right-hand-side may 
depend only on the independent variable t. As an example, consider a mass falling 
under the influence of constant gravity, such as approximately found on the Earth’s 
surface. Newton’s law, F = ma, results in the equation 


d2x 
ape — 8: 


where x is the height of the object above the ground, m is the mass of the object, and 
g = 9.8 meter / sec? is the constant gravitational acceleration. As Galileo suggested, 
the mass cancels from the equation, and 


a 
a * 


Here, the right-hand-side of the ode is a constant. The first integration, obtained by 
antidifferentiation, yields 


dx 
gp Oe hi 


with A the first constant of integration; and the second integration yields 
1 
x=B+At— =et’, 
2 
with B the second constant of integration. The two constants of integration A and 


B can then be determined from the initial conditions. If we know that the initial 
height of the mass is xg, and the initial velocity is vg, then the initial conditions are 


dx 
x0) =, (0) = 0. 
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1.1. THE SIMPLEST TYPE OF DIFFERENTIAL EQUATION 


Substitution of these initial conditions into the equations for dx/dt and x allows us 
to solve for A and B. The unique solution that satisfies both the ode and the initial 
conditions is given by 


x(t) = x9 + vot — 58 (1.1) 


For example, suppose we drop a ball off the top of a 50 meter building. How long 
will it take the ball to hit the ground? This question requires solution of (1.1) for 
the time T it takes for x(T) = 0, given x9 = 50 meter and vp = 0. Solving for T, 
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Chapter 2 


First-order differential 
equations 


Reference: Boyce and DiPrima, Chapter 2 


The general first-order differential equation for the function y = y(x) is written as 


4Y _ f(x,y), (2.1) 


where f(x,y) can be any function of the independent variable x and the dependent 
variable y. We first show how to determine a numerical solution of this equa- 
tion, and then learn techniques for solving analytically some special forms of (2.1), 
namely, separable and linear first-order equations. 


2.1 The Euler method 


View tutorial on YouTube 


Although it is not always possible to find an analytical solution of (2.1) for y = 
y(x), it is always possible to determine a unique numerical solution given an initial 
value y(xo) = yo, and provided f(x,y) is a well-behaved function. The differential 
equation (2.1) gives us the slope f(x9,yo) of the tangent line to the solution curve 
y = y(x) at the point (x9,yo). With a small step size Ax = x1 — xo, the initial 
condition (x9, Yo) can be marched forward to (x1,y1) along the tangent line using 
Euler’s method (see Fig. 2.1) 


¥1 = Yo + Axf (Xo, Yo). 


This solution (x1,¥1) then becomes the new initial condition and is marched for- 
ward to (x2, ¥2) along a newly determined tangent line with slope given by f (x1, 1). 
For small enough Ax, the numerical solution converges to the exact solution. 
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Yq = Vot AX F(X 5¥q) 


yy 


slope = f(X5:Y9) 


Figure 2.1: The differential equation dy/dx = f(x,y), y(xo) = yo, is integrated to x = x1 
using the Euler method y; = yo + Axf (xo, Yo), with Ax = x1 — x0. 


2.2 Separable equations 
View tutorial on YouTube 


A first-order ode is separable if it can be written in the form 


ety) = fl), y(20) = 40 2.2) 


where the function g(y) is independent of x and f(x) is independent of y. Integra- 
tion from xg to x results in 


[ suey! ax = a F(x)dx. 


The integral on the left can be transformed by substituting u = y(x), du = y'(x)dx, 
and changing the lower and upper limits of integration to y(x9) = yo and y(x) = y. 
Therefore, 


[sau = [faa 


and since uv is a dummy variable of integration, we can write this in the equivalent 
form 


[ “s(v)dy = f “flax. (2.3) 


A simpler procedure that also yields (2.3) is to treat dy/dx in (2.2) like a fraction. 
Multiplying (2.2) by dx results in 


g(y)dy = f(x)dx, 


which is a separated equation with all the dependent variables on the left-side, and 
all the independent variables on the right-side. Equation (2.3) then results directly 
upon integration. 
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Example: Solve ay + 5y = 3, with y(0) =2. 


We first manipulate the differential equation to the form 


dy 
ae =5(0- y), (2.4) 


and then treat dy/dx as if it was a fraction to separate variables: 


We integrate the right-side from the initial condition x = 0 to x and the left-side 
from the initial condition y(0) = 2 to y. Accordingly, 


. aT - -5/ dx. (2.5) 


The integrals in (2.5) need to be done. Note that y(x) < 3 for finite x or the integral 
on the left-side diverges. Therefore, 3 — y > 0 and integration yields 


Since this is our first nontrivial analytical solution, it is prudent to check our result. 
We do this by differentiating our solution: 


dy _1 1, 
dx 2 

1 

5 (3-9); 


and checking the initial conditions, y(0) = 3—e? = 2. Therefore, our solution 
satisfies both the original ode and the initial condition. 


Example: Solve ay + 5y = 3, with y(0) =4. 


This is the identical differential equation as before, but with different initial condi- 
tions. We will jump directly to the integration step: 


y aT _ 5 [ax 
43-4 2 
Now y(x) > 3, so that y—3 > 0 and an yields 
1 
—In(y—3)]q = 5]. 


In(y—3) = Es 
y—-3=0e°2%, 
y=3+e72%, 
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dy/dx + y/2 = 3/2 


Figure 2.2: Solution of the following ode: ay + 5y=3. 


The solution curves for a range of initial conditions is presented in Fig. 2.2. All 
solutions have a horizontal asymptote at y = 3 at which dy/dx = 0. For y(0) = yo, 
the general solution can be shown to be y(x) = 3 + (yo — 3) exp(—x/2). 

: dy _ 2cos2 
Example: Solve % = soy 
the solution exist? (ii) For what value of x > 0 is y(x) maximum? 


, with y(0) = —1. (i) For what values of x > 0 does 


Notice that the derivative of y diverges when y = —3/2, and that this may cause 
some problems with a solution. 
We solve the ode by separating variables and integrating from initial conditions: 


(3 + 2y)dy = 2cos2x dx 
[6+ 2ndy =2 [° cos2xax 
3yt+y"]", =sin2x]) 

y? + 3y+2-—sin2x =0 
i s[-3+ V1+4sin 2x]. 


Solving the quadratic equation for y has introduced a spurious solution that does 
not satisfy the initial conditions. We test: 


ys(0) = 5-34 1]={ 3 


Only the + root satisfies the initial condition, so that the unique solution to the ode 
and initial condition is 


pes 5[-3+ V1 + 4sin 2x]. (2.6) 


To determine (i) the values of x > 0 for which the solution exists, we require 


1+4sin2x > 0, 
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or j 
sin2x > — rt (2.7) 


Notice that at x = 0, we have sin2x = 0; at x = 71/4, we have sin2x = 1; at 
x = 7/2, we have sin2x = 0; and at x = 37/4, we have sin2x = —1 We therefore 
need to determine the value of x such that sin2x = —1/4, with x in the range 
m/2< x < 37/4. The solution to the ode will then exist for all x between zero and 
this value. 

To solve sin2x = —1/4 for x in the interval 7/2 < x < 37/4, one needs to 
recall the definition of arcsin, or sin~!, as found on a typical scientific calculator. 
The inverse of the function 


f(x) =sinx, -—7m/2<x< 7/2 


is denoted by arcsin. The first solution with x > 0 of the equation sin2x = —1/4 
places 2x in the interval (7,37/2), so to invert this equation using the arcsine 
we need to apply the identity sin(7— x) = sinx, and rewrite sin2x = —1/4 as 
sin (7 — 2x) = —1/4. The solution of this equation may then be found by taking 
the arcsine, and is 

m — 2x = arcsin (—1/4), 


or 


1 fool 
x= 3 (7+ arcsin +) . 
Therefore the solution exists for 0 < x < (m+ arcsin (1/4)) /2 = 1.6971..., where 
we have used a calculator value (computing in radians) to find arcsin(0.25) = 
0.2527.... At the value (x,y) = (1.6971...,—3/2), the solution curve ends and 
dy/dx becomes infinite. 

To determine (ii) the value of x at which y = y(x) is maximum, we examine 
(2.6) directly. The value of y will be maximum when sin2x takes its maximum 
value over the interval where the solution exists. This will be when 2x = 7/2, or 
x = 7/4 =0.7854.... 

The graph of y = y(x) is shown in Fig. 2.3. 


2.3 Linear equations 
View tutorial on YouTube 


The linear first-order differential equation (linear in y and its derivative) can be 
written in the form 


dy = 
dx + POX)Y = 8(2), (2.8) 


with the initial condition y(xo) = yo. Linear first-order equations can be integrated 
using an integrating factor p(x). We multiply (2.8) by p(x), 


wa) | + play] = wg, 2.) 


and try to determine ji(x) so that 
dy _d 
wx) | + playy| = Zln(aoa. (2.10) 
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(3+2y) dy/dx = 2 cos 2x, y(0) = —1 


0.5 1 1.5 
X 


Figure 2.3: Solution of the following ode: (3 + 2y)y’ = 2cos 2x, y(0) = —1. 


Equation (2.9) then becomes 


d 
ge HOY) = HX) 8 (x). (2.11) 
Equation (2.11) is easily integrated using pu(xo) = Wo and y(x9) = yo: 


x 


u(x)¥ — Hoyo = | : u(x)g(x)dx, 


or 
1 x 
= —_ + : x)e(x ax) ; 2.12 
Y= aeay (novo +f nldste) @.12) 
It remains to determine p(x) from (2.10). Differentiating and expanding (2.10) 
yields 
dy | — ah 4 WY. 
dx | PRY Gy ¥ 7 Bay? 
and upon simplifying, 
d 
oF _ py. (2.13) 


Equation (2.13) is separable and can be integrated: 


Hv du f 
SS x)dx, 
i H x0 a 


H x 
In— = / p(x)dx, 
Ho Xo 
x 
L(x) = Mo exp (/ pls)dx). 

xo 
Notice that since y is proportional to 9, the initial value fg cancels out of (2.12). It 
is therefore customary to simply assign fig = 1. The solution to (2.8) satisfying the 
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initial condition y(x9) = yo is then commonly written as 


Y= ata (vot [ma s(arde), 


n(x) = exp ( [ p(a)dr) 


the integrating factor. This important result finds frequent use in applied mathe- 
matics. 


with 


Example: Solve ay +2y =e *, with y(0) = 3/4. 


Note that this equation is not separable. With p(x) = 2 and g(x) = e*, we have 


u(x) = exp ([ 2a) 


= e*, 
and 
x. 
poe (24 [eterar 
0 
Fe & 
=e s+ ff dx) 
0 
i+-0) 


| 
N 
R 


II lI 

fan) ey 

N 

R 
NN Nom 


ray 
B 

rv | 
Ble 
Bigg 


Example: Solve ay — 2xy = x, with y(0) = 0. 


This equation is separable, and we solve it in two ways. First, using an integrating 
factor with p(x) = —2x and g(x) = x: 


u(x) = exp (-2 ye vd.) 


and 2 
y= aa: xe™ dx. 
0 
The integral can be done by substitution with u = x7, du = 2xdx: 


2 
4 P 4 
ih xe dx = | e “du 
0 0 


eek —u a 
=e | 


30-2") 
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Therefore, 


y dy _. f* 
coreg) xdx, 


zim (14+ 2y) = ae 


1+2y= ev, 


HE) 


The results from the two different solution methods are the same, and the choice of 
method is a personal preference. 


2.4 Applications 


2.4.1 Compound interest 


View tutorial on YouTube 


The equation for the growth of an investment with continuous compounding of 
interest is a first-order differential equation. Let S(t) be the value of the investment 
at time t, and let r be the annual interest rate compounded after every time interval 
At. We can also include deposits (or withdrawals). Let k be the annual deposit 
amount, and suppose that an installment is deposited after every time interval At. 
The value of the investment at the time f + At is then given by 


S(t + At) = S(t) + (rA#)S(t) + kAt, (2.14) 


where at the end of the time interval At, rAtS(t) is the amount of interest credited 
and kAt is the amount of money deposited (k > 0) or withdrawn (k < 0). Asa 
numerical example, if the account held $10,000 at time t, and r = 6% per year and 
k = $12,000 per year, say, and the compounding and deposit period is At = 1 month 
= 1/12 year, then the interest awarded after one month is rAtS = (0.06/12) x 
$10,000 = $50, and the amount deposited is kAt = $1000. 

Rearranging the terms of (2.14) to exhibit what will soon become a derivative, 
we have 

S(t + At) — S(t) 
At 
The equation for continuous compounding of interest and continuous deposits is 
obtained by taking the limit At — 0. The resulting differential equation is 
ds 


a ars tk, (2.15) 


= rS(t) +k. 
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which can solved with the initial condition S(0) = Sg, where So is the initial capital. 
We can solve either by separating variables or by using an integrating factor; I solve 
here by separating variables. Integrating from t = 0 to a final time t, 


S t 
/ dS z | dt, 
So tS +k 0 
1 rS+k 
oo (= = ae 
rS +k = (rSo + ke", 
ge: rSoe™ = = k 


S = Soe? + mon (l1-e"), (2.16) 


where the first term on the right-hand side of (2.16) comes from the initial invested 
capital, and the second term comes from the deposits (or withdrawals). Evidently, 
compounding results in the exponential growth of an investment. 

As a practical example, we can analyze a simple retirement plan. It is easiest to 
assume that all amounts and returns are in real dollars (adjusted for inflation). Sup- 
pose a 25 year-old plans to set aside a fixed amount every year of his/her working 
life, invests at a real return of 6%, and retires at age 65. How much must he/she 
invest each year to have HK$8,000,000 at retirement? (Note: US$1 ~ HK $8.) We 
need to solve (2.16) for k using t = 40 years, S(t) = $8,000,000, Sy = 0, and r = 0.06 
per year. We have 


rS(t) 
er ee hy 
k= 0.06 x 8,000,000 
_ e0.06x40 __ 4’ 


= $47,889 year o 


To have saved approximately one million US$ at retirement, the worker would need 
to save about HK$50,000 per year over his/her working life. Note that the amount 
saved over the worker’s life is approximately 40 x $50,000 = $2,000,000, while the 
amount earned on the investment (at the assumed 6% real return) is approximately 
$8,000,000 — $2,000,000 = $6,000,000. The amount earned from the investment is 
about 3x the amount saved, even with the modest real return of 6%. Sound invest- 
ment planning is well worth the effort. 


2.4.2 Chemical reactions 


Suppose that two chemicals A and B react to form a product C, which we write as 


ADRAC, 


where k is called the rate constant of the reaction. For simplicity, we will use the 
same symbol C, say, to refer to both the chemical C and its concentration. The law 
of mass action says that dC /dt is proportional to the product of the concentrations 
A and B, with proportionality constant k; that is, 


dC 
p= KAB. (2.17) 
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Similarly, the law of mass action enables us to write equations for the time-derivatives 
of the reactant concentrations A and B: 


— =-kAB, — =-—kAB. (2.18) 
The ode given by (2.17) can be solved analytically using conservation laws. We 


assume that Ag and Bo are the initial concentrations of the reactants, and that no 
product is initially present. From (2.17) and (2.18), 


“(A+C)=0 Ss Acs Ay 
d 
qB+C) =0 = B+C=Bpo. 


Using these conservation laws, (2.17) becomes 


- =k(Ap—C)(Bo-C), C(0) =0, 


which is a nonlinear equation that may be integrated by separating variables. Sep- 
arating and integrating, we obtain 


C dC t 
i amram 7th 
— kt. (2.19) 


The remaining integral can be done using the method of partial fractions. We write 


1 a b 


(Ao — C)(Bo — C) ~ Ag —C eae od (2.20) 


The cover-up method is the simplest method to determine the unknown coefficients 
a and b. To determine a, we multiply both sides of (2.20) by Ag — C and set C = Ag 
to find 


— 1 
~ By— Ao 
Similarly, to determine b, we multiply both sides of (2.20) by By — C and set C = Bo 
to find 1 
b= , 
Ao — Bo 


Therefore, 


1 _ 1 1 1 
(Ap —C)(Bo -—C)  Bo— Ao Ag—-C Bo—C/’ 
and the remaining integral of (2.19) becomes (using C < Apo, Bo) 


[ dC 1 ( Cc dC Cc dC ) 
0 (Ap — C)(Bo — C) Bo — Ao 0 Ag—C 0 Bo —C 
1 Ao-C)\ , ne 
ee in( Ao ) rin ( Bo )) 


1 Ao(Bo — C) 
By — Ao ™ (a = a 
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Using this integral in (2.19), multiplying by (Bp — Ag) and exponentiating, we obtain 


Ao(Bo — C) 


Bo—Ao)kt_ 
Bo(Ao — C) 


=p 


Solving for C, we finally obtain 


e(Bo—Ao)kt = 1 


Boe(Bo-Ao)kt — Ag’ 


C(t) = ApBo 


which appears to be a complicated expression, but has the simple limits 


lim C(t) = 


t—-yoo 


Bo, if Bo < Ao 
= min(Ao, Bo). 


is if Ay < Bo, 


As one would expect, the reaction stops after one of the reactants is depleted; and 
the final concentration of product is equal to the initial concentration of the depleted 
reactant. 


2.4.3 Terminal velocity 


View tutorial on YouTube 


Using Newton’s law, we model a mass m free falling under gravity but with air 
resistance. For simplicity, we assume that the force of air resistance is proportional 
to the speed of the mass and opposes the direction of motion. We define the x-axis to 
point in the upward direction, opposite the force of gravity. Near the surface of the 
Earth, the force of gravity is approximately constant and is given by —mg, with g = 
9.8 m/s? the usual gravitational acceleration. The force of air resistance is modeled 
by —kv, where v is the vertical velocity of the mass and k is a positive constant 
that depends on the size and shape of the falling mass. When the mass is falling, 
v < 0 and the force of air resistance is positive, pointing upward and opposing the 
motion. The total force on the mass is therefore given by F = —mg —kv. With 
F = ma and a = dv/dt, we obtain the differential equation 


m— = —mg —ko. (2.21) 


The terminal velocity vc of the mass is defined as the asymptotic velocity after air 
resistance balances the gravitational force. When the mass is at terminal velocity, 
dv/dt = 0 so that 


js. (2.22) 


The approach to the terminal velocity of a mass initially at rest is obtained by 
solving (2.21) with initial condition v(0) = 0. The equation is both linear and 
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separable, and I solve by separating variables: 


v du t 
mf ——— =- | dt, 
0 mg +kv 0 


ma (™*) — 
k mg 


14. =e kt/m 


mg 
v= 3 (1 — eH/m) . 


Therefore, v = Uco (1 ae ", and v approaches Ue as the exponential term de- 
cays to zero. 

As an example, a skydiver of mass m = 100kg with his parachute closed may 
have a terminal velocity of 200 km/hr. With 


g = (9.8m/s”)(10-3 km/m)(60s/min)?(60 min/hr)* = 127,008 km/hr’, 


one obtains from (2.22), k = 63,504k¢/hr. One-half of the terminal velocity for free- 
fall (100 km/hr) is therefore attained when (1 — g Him) = 1/2, ort = mln2/k & 
4sec. Approximately 95% of the terminal velocity (190km/hr ) is attained after 
17sec. 


2.4.4 Escape velocity 


View tutorial on YouTube 


An interesting physical problem is to find the smallest initial velocity for a mass 
on the Earth’s surface to escape from the Earth’s gravitational field, the so-called 
escape velocity. Newton’s law of universal gravitation asserts that the gravitational 
force between two massive bodies is proportional to the product of the two masses 
and inversely proportional to the square of the distance between them. For a mass 
m a position x above the surface of the Earth, the force on the mass is given by 


po Men 
(Rox)? 

where M and R are the mass and radius of the Earth and G is the gravitational 

constant. The minus sign means the force on the mass m points in the direction 

of decreasing x. The approximately constant acceleration g on the Earth’s surface 

corresponds to the absolute value of F/m when x = 0: 


_ GM 
sale 


and ¢ = 9.8 m/ s. Ignoring air resistance, Newton’s law F = ma for the mass m is 
thus given by 


d?x GM 
d2 ~~ (R+xp 
a & 
=-ay pe (2.23) 
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where the radius of the Earth is known to be R + 6350 km. 

A useful trick allows us to solve this second-order differential equation as a 
first-order equation. First, note that d?x/dt? = dv/dt. If we write v(t) = v(x(t))— 
considering the velocity of the mass m to be a function of its distance above the 
Earth—we have using the chain rule 


dv dudx 
dt dx dt 
_ dv 
=U he 


where we have used v = dx/dt. Therefore, (2.23) becomes the first-order ode 


dv g 


ax (1+x/R)?’ 


which may be solved assuming an initial velocity v(x = 0) = vp when the mass is 
shot vertically from the Earth’s surface. Separating variables and integrating, we 


obtain 
[ q [ dx 
v= ‘ 
. 8 fo (1+x/R)? 


The left integral is 3(e? — vj), and the right integral can be performed using the 
substitution vu =1+x«/R, du =dx/R: 


oe dx e 14+x/R gy 
I (+x/R)2 [ w 


x) 1+x/R 


Therefore, 


a) 2,  &Rx 
%) = x +R’ 
which when multiplied by m is an expression of the conservation of energy (the 
change of the kinetic energy of the mass is equal to the change in the potential 
energy). Solving for 0”, 
2.2 2gRx 
Uv" = 0 eR 
The escape velocity is defined as the minimum initial velocity vp such that the 
mass can escape to infinity. Therefore, 19 = Vescape When v + 0 as x — oo. Taking 
this limit, we have 


. 2@Rx 
a = tim, x+R 
= 2gR. 


VU, 


With R 6350 km and g = 127008 km/ hr’, we determine Vescape = \/2¢R ~ 40000 
km/hr. In comparison, the muzzle velocity of a modern high-performance rifle is 
4300 km/hr, almost an order of magnitude too slow for a bullet, shot into the sky, 
to escape the Earth’s gravity. 
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Figure 2.4: RC circuit diagram. 


2.4.5 RC circuit 


View tutorial on YouTube 


Consider a resister R and a capacitor C connected in series as shown in Fig. 2.4. 
A battery providing an electromotive force, or emf €, connects to this circuit by a 
switch. Initially, there is no charge on the capacitor. When the switch is thrown to 
a, the battery connects and the capacitor charges. When the switch is thrown to b, 
the battery disconnects and the capacitor discharges, with energy dissipated in the 
resister. Here, we determine the voltage drop across the capacitor during charging 
and discharging. 
The equations for the voltage drops across a capacitor and a resister are given 
by 
Vo = q/C, Vr =ik, (2.24) 


where C is the capacitance and R is the resistance. The charge q and the current i 
are related by 

. dg 

= 2.2 
i ai (2.25) 
Kirchhoff’s voltage law states that the emf € in any closed loop is equal to the 

sum of the voltage drops in that loop. Applying Kirchhoff’s voltage law when the 
switch is thrown to a results in 


Ve tVco=E. (2.26) 
Using (2.24) and (2.25), the voltage drop across the resister can be written in terms 
of the voltage drop across the capacitor as 
dVc 
Ve = RC—= 
R C ae 


and (2.26) can be rewritten to yield the linear first-order differential equation for Vc 
given by 
dVc 
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with initial condition Vc(0) = 0. 
The integrating factor for this equation is 


p(t) = elf RC, 
and (2.27) integrates to 
t 
Ve(t) =e RC f(E/RC)et/ Rat, 
0 
with solution 
Ve(t) = € (1- grunt) 


The voltage starts at zero and rises exponentially to €, with characteristic time scale 
given by RC. 
When the switch is thrown to b, application of Kirchhoff’s voltage law results in 


Vet+Vce =0, 
with corresponding differential equation 


ae +Vc/RC =0. 


Here, we assume that the capacitance is initially fully charged so that V-(0) = E. 
The solution, then, during the discharge phase is given by 


Vo(t) = Ee */RC, 
The voltage starts at € and decays exponentially to zero, again with characteristic 
time scale given by RC. 
2.4.6 The logistic equation 


View tutorial on YouTube 


Let N(t) be the number of individuals in a population at time ft, and let b and d be 
the average per capita birth rate and death rate, respectively. In a short time At, the 
number of births in the population is bAtN, and the number of deaths is dAtN. An 
equation for N at time ¢ + At is then determined to be 


N(t+ At) = N(t) + bDAtN(t) — dAtN(t), 
which can be rearranged to 


N(t+ At) — N(t) 


= (b—d)N(t); 
0 (b— a)N(W) 
and as At > 0, and with r = b — d, we have 
dN 
—=TN. 
dt" 


This is the Malthusian growth model (Thomas Malthus, 1766-1834), and is the same 
equation as our compound interest model. 
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Under a Malthusian growth model, the population size grows exponentially like 
N(t) = Noe", 


where Np is the initial population size. However, when the population growth is 
constrained by limited resources, a heuristic modification to the Malthusian growth 
model results in the Verhulst equation, 


dN N 
= =IN (1 = z) ; (2.28) 


where K is called the carrying capacity of the environment. Making (2.28) dimen- 


sionless using t = rt and x = N/K leads to the logistic equation, 


= x(1—~x), 


where we may assume the initial condition x(0) = x9 > 0. Separating variables and 


integrating 
x x 
| dx | dt. 
xo *¥ (1 _ x) 0 


The integral on the left-hand-side can be done using the method of partial fractions: 


1 a b 


x(l—x) x 1-2’ 


and the cover-up method yields a = b = 1. Therefore, 


a dx =| dx [ dx 
HX). dey Fda 2) 


=e in 

xo 1—x9 
_ 4, X= %0) 
g(t ae) 
=i 


Solving for x, we first exponentiate both sides and then isolate x: 


x(1 _ xo) = et 
xo(1—x) 
x(1— x0) = xpe’ — xxge", 
x(1— x9 + x9e") = xe", 
x0 
a ; 2.2 
** xo (1 — xo) — 


We observe that for xp > 0, we have lim;_,.. x(T) = 1, corresponding to 
lim N(t) = K. 
aoe 


The population, therefore, grows in size until it reaches the carrying capacity of its 
environment. 
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Linear second-order differential 
equations with constant 
coefficients 


Reference: Boyce and DiPrima, Chapter 3 


The general linear second-order differential equation with independent variable t¢ 
and dependent variable x = x(t) is given by 


¥+ p(t)x + q(t)x = g(t), (3.1) 


where we have used the standard physics notation x = dx/dt and ¥ = a2 x /dt?. 
A unique solution of (3.1) requires initial values x(tg) = xo and X(f9) = uo. The 
equation with constant coefficients—on which we will devote considerable effort— 
assumes that p(t) and q(t) are constants, independent of time. The second-order 
linear ode is said to be homogeneous if g(t) = 0. 


3.1 The Euler method 


View tutorial on YouTube 


In general, (3.1) cannot be solved analytically, and we begin by deriving an algo- 
rithm for numerical solution. Consider the general second-order ode given by 


¥ = f(t,x,X). 


We can write this second-order ode as a pair of first-order odes by defining u = x, 
and writing the first-order system as 


x=4U, (3.2) 
u = f(t,x,u). (3.3) 
The first ode, (3.2), gives the slope of the tangent line to the curve x = x(t); the 
second ode, (3.3), gives the slope of the tangent line to the curve u = u(t). Beginning 
at the initial values (x,u) = (xo,uo) at the time t = to, we move along the tangent 
lines to determine x; = x(to + At) and u, = u(tp + At): 
xX, = Xo + Atuo, 
uy = ug + Atf (to, Xo, Uo). 


The values x; and wu, at the time ft; = fo + At are then used as new initial values 
to march the solution forward to time ty = t; + At. As long as f(t,x,u) is a well- 
behaved function, the numerical solution converges to the unique solution of the 
ode as At > 0. 
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3.2 The principle of superposition 
View tutorial on YouTube 


Consider the linear second-order homogeneous ode: 

¥+ p(t)x + q(t)x =0; (3.4) 
and suppose that x = X;(t) and x = X(t) are solutions to (3.4). We consider a 
linear combination of X; and X> by letting 


X(t) = cy X(t) + c2Xo(t), (3.5) 


with c; and cz constants. The principle of superposition states that x = X(t) is also a 
solution of (3.4). To prove this, we compute 
X+ px + qx = cy Xy + cr Xo +p (c1Xy + co Xp) +4q (cy Xy + c2X2) 
Cl (X { pX t qX1) + C2 (Xp t pX2 t qX2) 
=c,x0+c¢. x0 
— 0, 


since X; and X» were assumed to be solutions of (3.4). We have therefore shown 
that any linear combination of solutions to the homogeneous linear second-order 
ode is also a solution. 


3.3. The Wronskian 


View tutorial on YouTube 


Suppose that having determined that two solutions of (3.4) are x = X (ft) and 
x = X(t), we attempt to write the general solution to (3.4) as (3.5). We must then 
ask whether this general solution will be able to satisfy the two initial conditions 
given by 

x(fo) = xo, (to) = uo. (3.6) 


Applying these initial conditions to (3.5), we obtain 


c1X1 (to) + ¢2Xo(to) = Xo, 
c1 Xi (to) + c2X2(to) = uo, (3.7) 


which is observed to be a system of two linear equations for the two unknowns c1 
and cz. Solution of (3.7) by standard methods results in 


_ XoX2(to) — uoX2(to) gies ugX1 (to) — xoX1 (fo) 
1 WwW y 2. WwW , 


where W is called the Wronskian and is given by 
W = Xi (to) X2(to) — X1(to) X2(to). (3.8) 


Evidently, the Wronskian must not be equal to zero (W ¥ 0) for a solution to exist. 
For examples, the two solutions 


X\(t) =Asinwt, X(t) = Bsinwt, 


32 CHAPTER 3. SECOND-ORDER ODES, CONSTANT COEFFICIENTS 


3.4. HOMOGENEOUS ODES 


have a zero Wronskian at f = tg, as can be shown by computing 
W = (Asinwtg) (Bw cos wtg) — (Aw cos wig) (B sin wt) 
= 0; 
while the two solutions 
X1(t) =sinwt, X(t) = coswt, 
with w 4 0, have a nonzero Wronskian at t = fo, 


W = (sinwtg) (—w sinwtg) — (w cos wt) (cos wtg) 


=—-W. 


When the Wronskian is not equal to zero, we say that the two solutions X(t) and 
X(t) are linearly independent. The concept of linear independence is borrowed 
from linear algebra, and indeed, the set of all functions that satisfy (3.4) can be 
shown to form a two-dimensional vector space. 


3.4 Homogeneous linear second-order ode with con- 
stant coefficients 


View tutorial on YouTube 


We now study solutions of the homogeneous, constant coefficient ode, written as 
ax + bx +cx =0, (3.9) 


with a, b, and c constants. Such an equation arises for the charge on a capacitor 
in an unpowered RLC electrical circuit, or for the position of a freely-oscillating 
frictional mass on a spring, or for a damped pendulum. Our solution method finds 
two linearly independent solutions to (3.9), multiplies each of these solutions by a 
constant, and adds them. The two free constants can then be used to satisfy two 
given initial conditions. 

Because of the differential properties of the exponential function, a natural 
ansatz, or educated guess, for the form of the solution to (3.9) is x = e!, where 
r is a constant to be determined. Successive differentiation results in x = re’ and 
% =7’e", and substitution into (3.9) yields 


are" + bre’ + ce =0. (3.10) 


Our choice of exponential function is now rewarded by the explicit cancelation in 
(3.10) of e”'. The result is a quadratic equation for the unknown constant r: 


ar? + br+c=0. (3.11) 


Our ansatz has thus converted a differential equation into an algebraic equation. 
Equation (3.11) is called the characteristic equation of (3.9). Using the quadratic for- 
mula, the two solutions of the characteristic equation (3.11) are given by 


1 

re = a (- ry b2 — 4ac) ‘ 

There are three cases to consider: (1) if b2 — 4ac > 0, then the two roots are distinct 
and real; (2) if b? — 4ac < 0, then the two roots are distinct and complex conjugates 
of each other; (3) if b? — 4ac = 0, then the two roots are degenerate and there is only 
one real root. We will consider these three cases in turn. 
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3.4.1 Distinct real roots 


When r+ # r— are real roots, then the general solution to (3.9) can be written as a 
linear superposition of the two solutions e+! and e’-'; that is, 


x(t) = cye"t* + eye’ 


The unknown constants c; and c2 can then be determined by the given initial con- 
ditions x(f9) = x9 and x(t) = uo. We now present two examples. 


Example 1: Solve ¥ + 5x + 6x = 0 with x(0) = 2, x(0) = 3, and find the maximum 
value attained by x. 


View tutorial on YouTube 
We take as our ansatz x = e’ and obtain the characteristic equation 
r+5r+6=0, 


which factors to 
(r+3)(r+2) =0. 


The general solution to the ode is thus 
x(t) = ce 7 + ene, 
The solution for x obtained by differentiation is 
z(t) = —2cye~2 — 3ene—*. 


Use of the initial conditions then results in two equations for the two unknown 
constant c, and cp: 


cy +2 = 2, 
—2cy — 3c2 =3. 


Adding three times the first equation to the second equation yields cj = 9; and the 
first equation then yields cp = 2—c, = —7. Therefore, the unique solution that 
satisfies both the ode and the initial conditions is 


a(f)=9e 7 =7e 


= 9e~2# (1 _ se) ; 


Note that although both exponential terms decay in time, their sum increases ini- 
tially since x(0) > 0. The maximum value of x occurs at the time t, when x = 0, 
or 


tm = In(7/6). 
The maximum xX = x(tm) is then determined to be 


Xm = 108/49. 
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Example 2: Solve ¢— x = 0 with x(0) = xo, x(0) = uo. 


Again our ansatz is x = e™ and we obtain the characteristic equation 


r?-1=0, 


with solution r; = +1. Therefore, the general solution for x is 


x(t) = cye' + c.e', 


and the derivative satisfies 
t 


x(t) = cye’ — coe. 
Initial conditions are satisfied when 


Cy + C2 = Xo, 


C1 — Co = Uo. 
Adding and subtracting these equations, we determine 
1 1 
C1 = 5 (Xo + uo), c= 5 (Xo — uo), 


so that after rearranging terms 


The terms in parentheses are the usual definitions of the hyperbolic cosine and sine 
functions; that is, 
t -t t —t 
e+e : e—e 
, sinht = 
2 2 


Our solution can therefore be rewritten as 


cosht = 


x(t) = x9 cosht + up sinht. 


Note that the relationships between the trigonometric functions and the complex 
exponentials were given by 


eit 4 eit git Ht 
cost = ————.,, sint = —__,, 
2. 21 
so that 
cosht = cosit, sinht = —isinit. 


Also note that the hyperbolic trigonometric functions satisfy the differential equa- 
tions 


La sinh t = cosht, fe cosht = sinht, 
dt dt 


which are similar to the differential equations satisfied by the more commonly used 
trigonometric function, but is absent a minus sign for the derivative of cosht. 
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3.4.2 Distinct complex-conjugate roots 
View tutorial on YouTube 


We now consider a characteristic equation (3.11) with b* — 4ac < 0, so the roots 
occur as complex-conjugate pairs. With 


= Z V 4ac — b2, 


A=-— = 
da’ + 99 


the two roots of the characteristic equation are A + ip and A — ip. We have thus 
found the following two complex exponential solutions to the differential equation: 


Z,(t) = eM! pitt Z(t) = eMto—ipe 


Applying the principle of superposition, any linear combination of Z; and Zp» is 
also a solution to the second-order ode. 

Recall that if z = x + iy, then Rez = (z+ Z)/2 and Imz = (z —2Z)/2i. We can 
therefore form two different linear combinations of Z(t) and Z(t) that are real, 
namely X;(t) = Re Z,(t) and X2(t) = Im Z;(t). We have 


X(t) =e" cost, Xo(t) =e sin pt. 


Having found these two real solutions, X;(t) and X(t), we can then apply the 
principle of superposition a second time to determine the general solution for x(t): 


x(t) =e (Acos pt + Bsinypt). (3.12) 


It is best to memorize this result. The real part of the roots of the characteristic equa- 
tion goes into the exponential function; the imaginary part goes into the argument 
of cosine and sine. 


Example 1: Solve ¢+ x = 0 with x(0) = xp and x(0) = uo. 


View tutorial on YouTube 


The characteristic equation is 


with roots 


r+ =i. 


The general solution of the ode is therefore 
x(t) = Acost+ Bsint. 


The derivative is 
x(t) = —Asint + Bcost. 


Applying the initial conditions: 
x(0) =A=xo, x(0) = B= up; 
so that the final solution is 
x(t) = xg cost + ug sint. 
Recall that we wrote the analogous solution to the ode *—x = 0 as x(t) = 


xq cosh t + ug sinh t. 
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Example 2: Solve #+ *%+ x =0 with x(0) = 1 and x(0) =0. 


The characteristic equation is 


rP+rt+1=0, 
with roots 

, = 1413 

eee ae a 


The general solution of the ode is therefore 


x(t) = ew! (Aco V3 44 Bsin #) : 


The derivative is 


1 
x(t) = —5e 2! (Aco V3 44 Bsin %) 


Applying the initial conditions x(0) = 1 and x(0) = 0: 


A=l1, 

—1 J3 

—A+—B=0; 

5) + 5 0; 
or 

Ani, pa. 

3 

Therefore, 


3.4.3 Repeated roots 

View tutorial on YouTube 

Finally, we consider the characteristic equation, 
ar>+br+c= 0, 


with b? — 4ac = 0. The degenerate root is then given by 


(3.13) 


To satisfy two initial conditions, a second independent solution must be found with 
nonzero Wronskian, and apparently this second solution is not of the form of our 


ansatz x = exp (rf). 
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One method to determine this missing second solution is to try the ansatz 


x(t) = y(t)xi(t), (3.14) 


where y(f) is an unknown function that satisfies the differential equation obtained 
by substituting (3.14) into (3.9). This standard technique is called the reduction of 
order method and enables one to find a second solution of a homogeneous linear 
differential equation if one solution is known. If the original differential equation is 
of order n, the differential equation for y = y(t) reduces to an order one lower, that 
is,n—1. 

Here, however, we will determine this missing second solution through a lim- 
iting process. We start with the solution obtained for complex roots of the charac- 
teristic equation, and then arrive at the solution obtained for degenerate roots by 
taking the limit p — 0. 

Now, the general solution for complex roots was given by (3.12), and to properly 
limit this solution as  — 0 requires first satisfying the specific initial conditions 
x(0) = x9 and x(0) = wo. Solving for A and B, the general solution given by (3.12) 
becomes the specific solution 


x(t;y) =e (+ cos pt + ae sin ut) , 


Here, we have written x = x(t; j1) to show explicitly that x depends on pi. 
Taking the limit as » — 0, and using lim, yu! sin ut = t, we have 

lim x(t; =e (x9 + (ug — Axo)E). 

in (tH) (xo + (wo o)f) 
The second solution is observed to be a constant, U9 — Axg, times t times the first 
solution, e“!. Our general solution to the ode (3.9) when b? — 4ac = 0 can therefore 
be written in the form 

x(t) = (cy + cot)e"™, 

where r is the repeated root of the characteristic equation. The main result to be 
remembered is that for the case of repeated roots, the second solution is ¢t times the 
first solution. 


Example: Solve «+ 2x + x = 0 with x(0) = 1 and x(0) = 0. 
The characteristic equation is 
P+2r+1=(r+1)P 
— 0, 
which has a repeated root given by r = —1. Therefore, the general solution to the 
ode is 
x(t) = ce * + egte, 


with derivative 
t 


X(t) = —cye~' + coe * — cyte !. 
Applying the initial conditions, we have 
c= 1, 
—Cy +c, =0, 
so that cy) = cp = 1. Therefore, the solution is 


ssi ne. 
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3.5 Inhomogeneous linear second-order ode 


We now consider the general inhomogeneous linear second-order ode (3.1): 
¥+ p(t)x + q(t)x = g(t), (3.15) 


with initial conditions x(t9) = xo and x(t9) = up. There is a three-step solution 
method when the inhomogeneous term g(t) # 0. (i) Find the general solution of 
the homogeneous equation 


¥+ p(t)x+q(t)x =0. (3.16) 
Let us denote the homogeneous solution by 
Xn (t) = c1X1(t) + c2Xa(t), 


where Xj and Xp are linearly independent solutions of (3.16), and cy and c2 are as 
yet undetermined constants. (ii) Find any particular solution x, of the inhomoge- 
neous equation (3.15). A particular solution is readily found when p(t) and q(t) are 
constants, and when ¢(f) is a combination of polynomials, exponentials, sines and 
cosines. (iii) Write the general solution of (3.15) as the sum of the homogeneous 
and particular solutions, 

x(t) = xn(t) + Xp(t), (3.17) 


and apply the initial conditions to determine the constants cy and cj. Note that 
because of the linearity of (3.15), 
7 1 y d? | | d | | | 
E+ pt + gx = so (%n + Xp) + Path + Xp) + (Xn + Xp) 
= (¥n + px + 9Xn) + (Fp + pXp + 9Xp) 
=0+8 
— 8, 
so that (3.17) solves (3.15), and the two free constants in x; can be used to satisfy 
the initial conditions. 


We will consider here only the constant coefficient case. We now illustrate the 
solution method by an example. 


Example: Solve ¢ — 3x — 4x = 3e”! with x(0) = 1 and x(0) = 0. 


View tutorial on YouTube 


First, we solve the homogeneous equation. The characteristic equation is 


r —3r—4= (r—4)(r+1) 
— 0, 
so that 
x,(t) = cre + cet. 


Second, we find a particular solution of the inhomogeneous equation. The form of 
the particular solution is chosen such that the exponential will cancel out of both 
sides of the ode. The ansatz we choose is 


x(t) = Ae”, (3.18) 
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where A is a yet undetermined coefficient. Upon substituting x into the ode, differ- 
entiating using the chain rule, and canceling the exponential, we obtain 


4A —6A-—4A =3, 


from which we determine A = —1/2. Obtaining a solution for A independent of ¢ 
justifies the ansatz (3.18). Third, we write the general solution to the ode as the sum 
of the homogeneous and particular solutions, and determine c; and cz that satisfy 
the initial conditions. We have 


1 
x(t) = cye#! + coe! — xe 


and taking the derivative, 


&(t) = 4e,e" — cge* — e*. 


Applying the initial conditions, 


or 


3 
cy + c2 = 57 
4cy —C2 = 1. 


This system of linear equations can be solved for cy by adding the equations to 
obtain cy = 1/2, after which cp = 1 can be determined from the first equation. 
Therefore, the solution for x(t) that satisfies both the ode and the initial conditions 
is given by 


1 1 
x(t) = se _ xf +e t 


= x (1 age + 20-5) : 


where we have grouped the terms in the solution to better display the asymptotic 
behavior for large ft. 

We now find particular solutions for some relatively simple inhomogeneous 
terms using this method of undetermined coefficients. 


Example: Find a particular solution of ¥ — 3x — 4x = 2sint. 
View tutorial on YouTube 
We show two methods for finding a particular solution. The first more direct 


method tries the ansatz 
x(t) = Acost+ Bsint, 


where the argument of cosine and sine must agree with the argument of sine in the 
inhomogeneous term. The cosine term is required because the derivative of sine is 
cosine. Upon substitution into the differential equation, we obtain 


(—Acost — Bsint) —3(—Asint + Bcost) —4(Acost+ Bsint) = 2sint, 
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or regrouping terms, 
— (5A +3B)cost+ (3A — 5B) sint = 2sint. 


This equation is valid for all t, and in particular for t = 0 and 7/2, for which the 
sine and cosine functions vanish. For these two values of t, we find 


5A+3B=0, 3A—5B=2; 


and solving for A and B, we obtain 


The particular solution is therefore given by 


1 
Xp = a5 (3cost —5sint). 


The second solution method makes use of the relation ec! = cost +isint to 


convert the sine inhomogeneous term to an exponential function. We introduce the 
complex function z(t) by letting 


z(t) = x(t) + y(t), 


where x(f) and y(t) are real functions, and rewrite the differential equation in com- 
plex form. We can rewrite the equation in one of two ways. On the one hand, if we 
use sint = Re{—ie"}, then the differential equation is written as 


Z—32-—4z = —2ie''. (3.19) 


and by equating the corresponding real and imaginary parts of (3.19), this complex 
equation becomes the two real differential equations 


X—3x—4x=2sint, yj —3y—4y = —2cost. 


The solution we are looking for, then, is xp(t) = Re{zp(t)}. 
On the other hand, if we write sint = Im{e''}, then the complex differential 
equation becomes 
Z—32—4z =2¢!!, (3.20) 


which becomes the two real differential equations 
¥—3x—4x=2cost, i —3y—4y =2sint. 


Here, the solution we are looking for is xp»(t) = Im{zp(t)}. 
We will proceed here by solving (3.20). As we now have an exponential function 
as the inhomogeneous term, we can make the ansatz 


a(t) = Ce", 
where we now expect C to be a complex constant. Upon substitution into the ode 
(3.20) and using i* = —1: 
—C—3iC —4C =2; 


CHAPTER 3. SECOND-ORDER ODES, CONSTANT COEFFICIENTS 41 


3.6. INHOMOGENEOUS LINEAR FIRST-ORDER ODES REVISITED 


or solving for C: 


—2 
aa ET 
—2(5 — 3i) 
(5 + 3i)(5 — 3i) 
—10+ 61 
34 
—5+4+3i 


17 


Therefore, 


= im { Z>(-5+3i)(cost-+isin4)} 


3cost—5sint). 


a7 
Example: Find a particular solution of # +x — 2x = f°. 
View tutorial on YouTube 
The correct ansatz here is the polynomial 

x(t) = At* + BE+C. 
Upon substitution into the ode, we have 


2A+2At+ B—2At? —2Bt—-2C =?, 


or 
DAP L2(A— BEL OA+ Ba ace =P. 


Equating powers of t, 
-2A=1, 2(A—B)=0, 2A+B-—2C=0; 


and solving, 
1 1 3 
A= B=— =—-, 
2’ 2’ . 4 


The particular solution is therefore 


3.6 Inhomogeneous linear first-order odes revisited 
The linear first-order ode can be solved by use of an integrating factor. Here I show 
that odes having constant coefficients can be solved by our newly learned solution 


method. 
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Example: Solve x + 2x = e~' with x(0) = 3/4. 


Rather than using an integrating factor, we follow the three-step approach: (i) find 
the general homogeneous solution; (ii) find a particular solution; (iii) add them 
and satisfy initial conditions. Accordingly, we try the ansatz x;,(t) = e for the 
homogeneous ode x + 2x = 0 and find 


r+2=0, or r=-2. 
To find a particular solution, we try the ansatz x,(f) = Ae‘, and upon substitution 
—-A+2A=1, or A=1. 
Therefore, the general solution to the ode is 
x(i)=ce te, 
The single initial condition determines the unknown constant c: 


x(0)=F=ct1, 


so that c = —1/4. Hence, 


3.7 Resonance 
View tutorial on YouTube 


Resonance occurs when the frequency of the inhomogeneous term matches the 
frequency of the homogeneous solution. To illustrate resonance in its simplest em- 
bodiment, we consider the second-order linear inhomogeneous ode 


¥+upx = fcoswt, x(0) =x, x(0) = uo. (3.21) 


Our main goal is to determine what happens to the solution in the limit w — wy. 
The homogeneous equation has characteristic equation 


r’+ur=0, 


so that r- = +iwp. Therefore, 


xp(t) = cy cos wot + cz sin wot. (3.22) 
To find a particular solution, we note the absence of a first-derivative term, and 


simply try 
x(t) = Acoswt. 


Upon substitution into the ode, we obtain 
—wAt+ weA =f; 
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or 


Therefore, 


Our general solution is thus 


x(t) = cy cos wot + c2 sin wot + —-~—, cos wt, 
wo — 
with derivative 
X(t) = wo(cz cos wot — cy sin wot) — zi sin wt 
wy — w 
Initial conditions are satisfied when 
xo = cy + = 
we — Ww" 
Ug = C2W9, 
so that 
Cy = Xo Z Q= Ho 
we — Ww?" wo 


Therefore, the solution to the ode that satisfies the initial conditions is 


f uo .. i 
t)= i nd Wot + — Wot + —— wt 
x(t) (x - x | COs wot + a sin Wot + me 2 cos 


up. cos wt — cos wot 
= x9 cos wot + — sinwot t ft 7 of) 
wo wy — w 


, 


where we have grouped together terms proportional to the forcing amplitude f. 

Resonance occurs in the limit w — wo; that is, the frequency of the inhomoge- 
neous term (the external force) matches the frequency of the homogeneous solution 
(the free oscillation). By L’Hospital’s rule, the limit of the term proportional to f is 
found by differentiating with respect to w: 


lim f(coswt —coswot) _ ise —ftsinwt 
2542 —2 
W->Wo9 WwW. W W->Wo9 WwW 
i (3.23) 
_ ftsinwot 
i 2wW9o ; 


At resonance, the term proportional to the amplitude f of the inhomogeneous term 
increases linearly with t, resulting in larger-and-larger amplitudes of oscillation for 
x(t). In general, if the inhomogeneous term in the differential equation is a solution 
of the corresponding homogeneous differential equation, then the correct ansatz for 
the particular solution is a constant times the inhomogeneous term times t. 

To illustrate this same example further, we return to the original ode, now as- 
sumed to be exactly at resonance, 


¥+ wax = f cos uot, 


44 CHAPTER 3. SECOND-ORDER ODES, CONSTANT COEFFICIENTS 


3.7. RESONANCE 


and find a particular solution directly. The particular solution is the real part of the 
particular solution of 

Sans = feo, 
If we try Zp = Ceo, we obtain 0 = f, showing that the particular solution is not 
of this form. Because the inhomogeneous term is a solution of the homogeneous 
equation, we should take as our ansatz 


_ iwot 
Zp = Ate, 


We have ; 
fp = Ac0! (1+ iwot), Zp = Ae (2iwo = w9t) ; 


and upon substitution into the ode 


Zp + WwaZp = Aeot (2ico — w§t) + wp Atelo! 


= 2iwy Aci! 
= felwot, 
Therefore, 
ri 
A=, 
2iwg 
and 


ft iwot 
= Ref ——. 0 
ss loa } 


_ ftsin wot 
~ 2W9 


the same result as (3.23). 


Example: Find a particular solution of ¢— 3x —4x =5e'. 


View tutorial on YouTube 


If we naively try the ansatz 
x= Aet, 


and substitute this into the inhomogeneous differential equation, we obtain 
A+3A-—4A=5, 


or 0 = 5, which is clearly nonsense. Our ansatz therefore fails to find a solution. The 
cause of this failure is that the corresponding homogeneous equation has solution 
i= cye"*# + coe, 


so that the inhomogeneous term 5e~‘ is one of the solutions of the homogeneous 
equation. To find a particular solution, we should therefore take as our ansatz 


x = Ate, 
with first- and second-derivatives given by 


a=Ae*(1—f, #= Ae '(—2++#). 
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Substitution into the differential equation yields 
Ae~'(—24+ t) —3Ae*(1 — t) —4Ate™' = 5e?. 
The terms containing t cancel out of this equation, resulting in -5A = 5, or A = —1. 
Therefore, the particular solution is 
Xp = —te!. 


3.8 Applications 


View Nondimensionalization on YouTube 


3.8.1 RLC circuit 
View RLC circuit on YouTube 


E(t) 


Figure 3.1: RLC circuit diagram. 


Consider a resister R, an inductor L, and a capacitor C connected in series as shown 
in Fig. 3.1. An AC generator provides a time-varying electromotive force (emf), 
E(t), to the circuit. Here, we determine the differential equation satisfied by the 
charge on the capacitor. 

The constitutive equations for the voltage drops across a capacitor, a resister, 
and an inductor are given by 


di 
Vo=a/e- Vek. Ve = <L, (3.24) 
where C is the capacitance, R is the resistance, and L is the inductance. The charge 
q and the current i are related by 
. dg 
el 2. 
i ai (3.25) 
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Kirchhoff’s voltage law states that the emf € applied to any closed loop is equal 
to the sum of the voltage drops in that loop. Applying Kirchhoff’s voltage law to 
the RLC ciruit results in 


Vi + Vrat+Ve = E(t); (3.26) 
or using (3.24) and (3.25), 


d’q dq 1 
L +R t = €(t). 
ae ae ete 
The equation for the RLC circuit is a second-order linear inhomogeneous differen- 

tial equation with constant coefficients. 
The AC voltage can be written as E(t) = Ey coswt, and the governing equation 
for q = q(t) can be written as 


Rdq | 1 _ & 
ge La tC! OL 


cos wt. (3.27) 


Nondimensionalization of this equation will be shown to reduce the number of free 
parameters. 

To construct dimensionless variables, we first define the natural frequency of 
oscillation of a system to be the frequency of oscillation in the absence of any driving 
or damping forces. The iconic example is the simple harmonic oscillator, with 
equation given by 

K+ wax =, 


and general solution given by x(t) = Acoswot + Bsinwot. Here, the natural fre- 
quency of oscillation is wo, and the period of oscillation is T = 27t/wy. For the 
RLC circuit, the natural frequency of oscillation is found from the coefficient of the 
q term, and is given by 
- 1 
=>: 
VLC 

Making use of wo, with units of one over time, we can define a dimensionless time 
t and a dimensionless charge Q by 


The resulting dimensionless equation for the RLC circuit is then given by 


dQ _ dQ _ 
gn ea +Q=cos BT, (3.28) 


where a and f are dimensionless parameters given by 


Notice that the original five parameters in (3.27), namely R, L, C, Ey) and w, have 
been reduced to the two dimensionless parameters a and f. We will return later to 
solve (3.28) after visiting two more applications that will be shown to be governed 
by the same dimensionless equation. 
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Figure 3.2: Mass-spring system (top view). 


3.8.2. Mass on a spring 


View Mass on a Spring on YouTube 


Consider a mass lying on a flat surface connected by a spring to a wall, with top 
view shown in Fig. 3.2. The spring force is modeled by Hooke’s law, F; = —kx, and 
sliding friction is modeled as Fr = —cdx/dt. An external force is applied and is 
assumed to be sinusoidal, with F. = Fy coswt. Newton’s equation, F = ma, results 
in 


dx dx 
m = —kx—-c + Fy cos wt. 


Rearranging terms, we obtain 


d’x  cdx | k_ fo - 
d2 ‘mdt' mm" 


Here, the natural frequency of oscillation is given by 


k 
wo — m’ 


and we define a dimensionless time t and a dimensionless position X by 


=> ta—=—+X=cos ft, (3.29) 


Notice that even though the physical problem is different, and the dimensionless 
variables are defined differently, the resulting dimensionless equation (3.29) for the 
mass-spring system is the same as that for the RLC circuit (3.28). 
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Figure 3.3: The pendulum. 


3.8.3 Pendulum 


View Pendulum on YouTube 


Here, we consider a mass that is attached to a massless rigid rod and is con- 
strained to move along an arc of a circle centered at the pivot point (see Fig. 3.3). 
Suppose ! is the fixed length of the connecting rod, and @ is the angle it makes with 
the vertical. 

We can apply Newton’s equation, F = ma, to the mass with origin at the bottom 
and axis along the arc with positive direction to the right. The position s of the mass 
along the arc is given by s = /6. The relevant gravitational force on the pendulum 
is the component along the arc, and from Fig. 3.3 is observed to be 


Fo = —mgsiné. 


We model friction to be proportional to the velocity of the pendulum along the arc, 
that is 
Fy = —c8 = —cld. 
With a sinusoidal external force, F, = Fy cos wt, Newton’s equation m& = Fy + Fr + 
F, results in 
mlé = —mg sin 6 — cl6 + Fy coswt. 
Rewriting, we have 
& 


, : Fi 
O+ aes © sing — — coswt. 
m ml 


At small amplitudes of oscillation, we can approximate sin @ ~ 6, and the natural 
frequency of oscillation wo of the mass is given by 


w= 8 
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Nondimensionalizing time as T = wot, the dimensionless pendulum equation be- 
comes 

a0 do, g= 

7 + oe + sin@ = ycos Bt, 


where a, 6, and are dimensionless parameters given by 


c Ww Fo 


“Kg = — => — — . 
mw’ p wo’ % mls 


The nonlinearity of the pendulum equation, with the term sin @, results in the addi- 
tional dimensionless parameter y. For small amplitude of oscillation, however, we 
can scale @ by 6 = 7@, and the small amplitude dimensionless equation becomes 


=> +a—=— +0 = cos Bt, (3.30) 


the same equation as (3.28) and (3.29). 


3.9 Damped resonance 
View Damped Resonance on YouTube 


We now solve the dimensionless equations given by (3.28), (3.29) and (3.30), written 
here as 


¥+ax+x=cos ft, (3.31) 


where the physical constraints of our three applications requires that a > 0. The 
homogeneous equation has characteristic equation 


r+ar+1=0, 
so that the solutions are 


fe ed 


2° 9 


When a? — 4 < 0, the motion of the unforced oscillator is said to be underdamped; 
when «2 — 4 > 0, overdamped; and when a2-4=0, critically damped. For all 
three types of damping, the roots of the characteristic equation satisfy Re(r+) < 0. 
Therefore, both linearly independent homogeneous solutions decay exponentially 
to zero, and the long-time asymptotic solution of (3.31) reduces to the non-decaying 
particular solution. Since the initial conditions are satisfied by the free constants 
multiplying the decaying homogeneous solutions, the long-time asymptotic solu- 
tion is independent of the initial conditions. 

If we are only interested in the long-time asymptotic solution of (3.31), we can 
proceed directly to the determination of a particular solution. As before, we con- 
sider the complex ode 


Z+az+z= eiBt 
with xp = Re(z,). With the ansatz z, = Ae'P', we have 
—p°A+iapA+A=1, 
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or 
1 
(—) + inp 


= (4 = a2 ER) (a pr) inp ) . 


To determine x», we utilize the polar form of a complex number. The complex 


number z = x + iy can be written in polar form as z = re'?, where r = \/x2 + y? 
and tan@ = y/x. We therefore write 


(3.32) 


(1 — B?) — in = re’, 


r= (1-6)? +0262, tang = oe 


Using the polar form, A in (3.32) becomes 


with 


= 1 oi 
a ( oe ’ 


and xp = Re(Ae’‘P') becomes 


i 
- ( (1— 6°)? + =) Re {e (B+9)} 


1 
= ( aa =) cos (Bi + ¢). 


We conclude with a couple of observations about (3.33). First, if the forcing 
frequency w is equal to the natural frequency w 9 of the undamped oscillator, then 

= land A = 1/ia, and xp = (1/a)sint. The oscillator position is observed 
to be 71/2 out of phase with the external force, or in other words, the velocity of 
the oscillator, not the position, is in phase with the force. Second, the value of 6 
that maximizes the amplitude of oscillation is the value of 6 that minimizes the 
denominator of (3.33). To determine 6, we thus minimize the function ¢(B") with 
respect to oe where 


(3.33) 


g(B°) =(1- By +07. 
Taking the derivative of g with respect to p’ and setting this to zero to determine 
Bm yields 
I-68) re =o, 


| a2 a2 
Bn = dig A 


the last approximation valid if a << 1 and the dimensionless damping coefficient 
is small. We can interpret this result by saying that small damping slightly lowers 
the “resonance” frequency of the undamped oscillator. 


or 
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Chapter 4 


The Laplace transform 


Reference: Boyce and DiPrima, Chapter 6 


The Laplace transform is most useful for solving linear, constant-coefficient 
ode’s when the inhomogeneous term or its derivative is discontinuous. Although 
ode’s with discontinuous inhomogeneous terms can also be solved by adopting al- 
ready learned methods, we will see that the Laplace transform technique provides 
a simpler, more elegant solution. 


4.1 Definitions and properties of the forward and in- 
verse Laplace transforms 


View tutorial on YouTube 


The main idea is to Laplace transform the constant-coefficient differential equation 
for x(t) into a simpler algebraic equation for the Laplace-transformed function X(s), 
solve this algebraic equation, and then transform X(s) back into x(t). The correct 
definition of the Laplace transform and the properties that this transform satisfies 
makes this solution method possible. 

An exponential ansatz is used in solving homogeneous constant-coefficient odes, 
and the exponential function correspondingly plays a key role in defining the Laplace 
transform. The Laplace transform of f(t), denoted by F(s) = L{f(t)}, is defined 
by the integral transform 


F(s) = [ ” et FE dt, (4.1) 


The improper integral given by (4.1) diverges if f(t) grows faster than e*! for large 
t. Accordingly, some restriction on the range of s may be required to guarantee 
convergence of (4.1), and we will assume without further elaboration that these 
restrictions are always satisfied. 

The Laplace transform is a linear transformation. We have 


Llerfilt) + eof} = [oe (exfilt) + cafslt))at 


= [ ef, (t)dt + cp - e’' f(t) dt 
= cL{filt)} + coL{folt)}- 


There is also a one-to-one correspondence between functions and their Laplace 
transforms. A table of Laplace transforms can therefore be constructed and used to 
find both Laplace and inverse Laplace transforms of commonly occurring functions. 
Such a table is shown in Table 4.1 (and this table will be distributed with the exams). 
In Table 4.1, 1 is a positive integer. Also, the cryptic entries for u;(t) and d(t —c) 
will be explained later in §4.3. 
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Le 


3. et! 

4. ¢" 

5. tet 

6. sin bt 

7. cos bt 

8. e” sin bt 
9. e* cos bt 
10. tsin bt 
11. tcos bt 


12. u-(t) 


14. d(t—c) 
15. x(t) 


16. ¥(t) 


f(t) =L£{F(s)} 


13. uc(t) f(t —c) 


F(s) = L{f(t)} 


F(s— a) 


e “F(s) 
e: 


sX(s) — x(0) 


s*X(s) — sx(0) — x(0) 


Table 4.1: Table of Laplace Transforms 
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The rows of Table 4.1 can be determined by a combination of direct integration 
and some tricks. We first compute directly the Laplace transform of e” f(t) (line 1) 


Lier f(t} = [eset f(tyat 


= [ee f(t 


= F(s—a). 


We also compute directly the Laplace transform of 1 (line 2): 


L{1}= i e “dt 


1 _ fore) 
er st 
§ 0 


1 


S 


Now, the Laplace transform of e” (line 3) may be found using these two results: 


Lhe} = Lie .1} 
1 


(4.2) 
sa 


The transform of t” (line 4) can be found by successive integration by parts. A more 
interesting method uses Taylor series for e” and 1/(s — a). We have 


a eo (at)” 

eqery =f a 
at (4.3) 
yet} 

Also, with s > a, 
i... 4 
s—a_ s(1— 4) 
1S sayn 
=) (4) 


lI 
ih 
Sls 
+/ 3 

An 


Using (4.2), and equating the coefficients of powers of a in (4.3) and (4.4), results in 
line 4: 


‘i n! 
Lit eg es 


The Laplace transform of t’e"* (line 5) can be found from line 1 and line 4: 


aca n} 
i ery cl 
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The Laplace transform of sin bt (line 6) may be found from the Laplace transform 
of e” (line 3) using a = ib: 


L{sin bt} = Im{L{el*}} 


1 
=im{—} 


Similarly, the Laplace transform of cos bt (line 7) is 


L{cos bt} = Re{L{ei} } 


~ Be 


The transform of e” sin bt (line 8) can be found from the transform of sin bt (line 6) 
and line 1: 


b 
tos = : 
L{e* sin bt} = (s—a2 +B’ 
and similarly for the transform of e” cos bt: 
at s—a 
bt} = ——_.—__.. 
L{e™ cos bt} Goan 


The Laplace transform of ¢ sin bt (line 10) can be found from the Laplace transform 
of te” (line 5 with n = 1) using a = ib: 


L{tsin bt} = Im{ L{ tel} } 


oe 
(s — ib)? 
(s + ib)? 
= I ns 
a { (s2 +22 
a 2bs 
~ (s2 af b2)2 . 
Similarly, the Laplace transform of t cos bt (line 11) is 


L{tcos bt} = Re{ L{tel}} 
(s + ib)? 
=Re lar Py 
52 — b2 
~ F4PP 


We now transform the inhomogeneous, second-order, constant-coefficient ode 
for x = x(t), 
ax + bx + cx = g(t), 


making use of the linearity of the Laplace transform: 
al{x} + bL{x} + cL{x} = L{g}. 
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To determine the Laplace transform of x(t) (line 15) in terms of the Laplace trans- 
form of x(t) and the initial conditions x(0) and x(0), we define X(s) = L{x(t)}, 


and integrate 
[oe) 
| e edt 
0 


by parts. We let 


Therefore, 
co : 00 wo 
| e xara xe | ior s | e*'xdt 
0 0 
= sX(s) — x(0), 
where assumed convergence of the Laplace transform requires 


: —st _ 
jim x(t)e * =0. 


Similarly, the Laplace transform of <(t) (line 16) is determined by integrating 


[oe) 
| e' ¥dt 
0 


by parts and using the just derived result for the first derivative. We let 
ie" dv = xdt 
du = —se~*'dt v=X, 
so that 
e dt = xe) +s [ e edt 
0 0 
—x(0) + s(sX(s) — x(0)) 
= s*X(s) — sx(0) — x(0), 


where similarly we assume lim;-00 X(t)e~! = 0. 


4.2 Solution of initial value problems 


We begin with a simple homogeneous ode and show that the Laplace transform 
method yields an identical result to our previously learned method. We then apply 
the Laplace transform method to solve an inhomogeneous equation. 


Example: Solve « — x — 2x = 0 with x(0) = 1 and x(0) = 0 by two different 
methods. 


View tutorial on YouTube 


The characteristic equation of the ode is determined from the ansatz x = e and is 
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The general solution of the ode is therefore 
x(t) = ce + coe. 


To satisfy the initial conditions, we must have 1 = c, +c2 and 0 = 2c, — co, requiring 
cq = 5 and c2 = 3. Therefore, the solution to the ode that satisfies the initial 


conditions is given by 
1 2 
a) = 3° ~ ae (4.5) 
We now solve this example using the Laplace transform. Taking the Laplace 


transform of both sides of the ode, using the linearity of the transform, and applying 
our result for the transform of the first and second derivatives, we find 


[s*X(s) — sx(0) — #(0)] — [sX(s) — x(0)] — [2X(s)] = 0, 


or 


(s — 1)x(0) + x(0) 

s*—s—2 ; 
Note that the denominator of the right-hand-side is just the quadratic from the 
characteristic equation of the homogeneous ode, and that this factor arises from the 
derivatives of the exponential term in the Laplace transform integral. 

Applying the initial conditions, we find 

s—1 
(s—2)(s +1) 
We have thus determined the Laplace transformed solution X(s) = L{x(t)}. We 
now need to compute the inverse Laplace transform x(t) = MX (s)}. 

However, direct inversion of (4.6) by searching Table 4.1 is not possible, but a 
partial fraction decomposition may be useful. In particular, we write 

s—1 #4 b 
S-VG4i)~ #2" 221" 

The cover-up method can be used to solve for a and b. We multiply both sides of 
(4.7) by s — 2 and put s = 2 to isolate a: 


X(s) = 


X(s) = 


(4.6) 


(4.7) 


_ eal 

ae q sa 
1 
=i! 

Similarly, we multiply both sides of (4.7) by s + 1 and put s = —1 to isolate b: 
s—1 

i s— ;| ae 
2 
=F 


Therefore, 
1 1 2 1 
X(s) = _- poi , 
(6) =3°5-9 +3 sa1 
and line 3 of Table 4.1 gives us the inverse transforms of each term separately to 
yield 


x(t) = 5a + at 
identical to (4.5). 
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Example: Solve ¥+ x = sin2t with x(0) = 2 and x(0) = 1 by Laplace transform 
methods. 


Taking the Laplace transform of both sides of the ode, we find 
s°X(s) —sx(0)— #(0) 4+ X(s) = C{sin 24} 
2 
~ s2 44’ 


where the Laplace transform of sin 2t made use of line 6 of Table 4.1. Substituting 
for x(0) and x(0) and solving for X(s), we obtain 

25-1 2 

~ s2+1 (s2+1)(s? +4)’ 


To determine the inverse Laplace transform from Table 4.1, we perform a partial 
fraction decomposition of the second term: 


2 as+b cs+d 


(s? + 1)(s2 + 4) = s2+1 s2 +40 (4.8) 


By inspection, we can observe that a = c = 0 and that d = —b. A quick calculation 
shows that 3b = 2, or b = 2/3. Therefore, 


2s+1 2/3 2/3 
X(s) = 
(s) se+1  s*+1  (s*+4) 
_ 2s 5/8 2/8 
st+1' s24+1  (s?44)’ 


From lines 6 and 7 of Table 4.1, we obtain the solution by taking inverse Laplace 
transforms of the three terms separately, where b = 1 in the first two terms, and 
b = 2 in the third term: 


1 
x(t) =2cost+ > sint — 3 sin 2t. 


4.3 Heaviside and Dirac delta functions 


The Laplace transform technique becomes truly useful when solving odes with 
discontinuous or impulsive inhomogeneous terms, these terms commonly modeled 
using Heaviside or Dirac delta functions. We will discuss these functions in turn, 
as well as their Laplace transforms. 


4.3.1 Heaviside function 


View tutorial on YouTube 


The Heaviside or unit step function (see Fig. 4.1) , denoted here by u,(t), is zero for 
t < c and is one for t > c; that is, 


0, t<c 
u-(t) = { 1, #oe (4.9) 


The precise value of u,(t) at the single point t = c shouldn’t matter. 
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x= u(t) 


t 


Figure 4.1: The Heaviside function. 


The Heaviside function can be viewed as the step-up function. The step-down 
function—one for t < c and zero for t > c—is defined as 


1, t<cG 
1—u;-(t) = { 0. te (4.10) 


The step-up, step-down function—zero for t < a, one for a < t < b, and zero for 
t > b—is defined as 


0, t<a@; 
ug(t)-—up(t)=4 1, a<t<b; (4.11) 
0, t>b. 


The Laplace transform of the Heaviside function is determined by integration: 


L{ue(t)} = i. eu, (t)dt 


[oe] 
— 1 edt 
c 
e 6 


== , 


S 


and is given in line 12 of Table 4.1. 


The Heaviside function can be used to represent a translation of a function f(t) 
a distance c in the positive t direction. We have 


0, t<c 


uc(t)f(t—c) = { f(t-c), t>c. 
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x=f(t) 


Figure 4.2: A linearly increasing function which turns into a constant function. 


The Laplace transform is 


Lluc(t)f(t—e)} = [oe stuclf(e— eda 
-[ve ~S F(t — c)dt 


= Se s(t! TeNe (t! (t") d t! 
= 7. eS! FE) dt! 
0 

=e “F(s), 
where we have changed variables to t' =t—c. The translation of f(t) a distance c in 
the positive t direction corresponds to the multiplication of F(s) by the exponential 
e—®. This result is shown in line 13 of Table 4.1. (Note that line 12 of Table 4.1 is a 
special case of line 13, where f(t) = 1.) 

Piecewise-defined inhomogeneous terms can be modeled using Heaviside func- 


tions. For example, consider the general case of a piecewise function defined on 
two intervals: 


MO={ AG, ase 
Using the Heaviside function u,, the function f(t) can be written in a single line as 
F(t) = fi) + (al) — falt))uc(#)- 
This example can be generalized to piecewise functions defined on multiple inter- 
ae a concrete example, suppose the inhomogeneous term is represented by a 


linearly increasing function, which then turns into a constant function, as sketched 


CHAPTER 4. THE LAPLACE TRANSFORM 61 


4.3. HEAVISIDE AND DIRAC DELTA FUNCTIONS 


in Fig. 4.2. Explicitly, 


t, ift<1; 


fi ={ 1 iff>1 


We can rewrite f(t) using the Heaviside function u(t): 


f(t) =t-m(t)(t-1); 


and we can take the Laplace transform of the function f(t) by writing 


F(s) = L{t} — L{uy(t) (t-1)}- 


The Laplace transform of the first term is found from line 4 of the table, and the 
Laplace transform of the second term is found from a combination of line 13 and 
line 4. We have 


4.3.2 Dirac delta function 


View tutorial on YouTube 


The Dirac delta function, denoted as 6(t), is defined by requiring that for any func- 


tion f(t), 
[ fOsat = F00). 


The usual view of the shifted Dirac delta function 6(t — c) is that it is zero every- 
where except at ¢ = c, where it is infinite, and the integral over the Dirac delta 
function is one. The Dirac delta function is technically not a function, but is what 
mathematicians call a distribution. Nevertheless, in most cases of practical interest, 
it can be treated like a function, where physical results are obtained following a 
final integration. 

There are many ways to represent the Dirac delta function as a limit of a well- 
defined function. For our purposes, the most useful representation makes use of 
the step-up, step-down function of (4.11) (see Fig. 4.3): 


1 
d(t-—c)= lim 5 (ue-el(t) — Uc+e(t)). 
Before taking the limit, the well-defined step-up, step-down function is zero except 
over a small interval of width 2e centered at tf = c, over which it takes the large 
value 1/2e. The integral of this function is one, independent of the value of e. 
The Laplace transform of the Dirac delta function is easily found by integration 
using the definition of the delta function: 


feeb a e-'§(t —c)dt 


— eS. 
This result is shown in line 14 of Table 4.1. 
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1/2 ¢ 


C-e Cte 


Figure 4.3: The Dirac delta function constructed using the step-up, step-down function. 


4.4 Discontinuous or impulsive inhomogeneous terms 


We now solve some more challenging ode’s with discontinuous or impulsive inho- 
mogeneous terms. 
Example: Solve 2% + % + 2x = us(t) — ug9(t), with x(0) = x(0) = 0 


The inhomogeneous term here is a step-up, step-down function that is unity over 
the interval (5,20) and zero elsewhere. Taking the Laplace transform of the ode 
using Table 4.1, 


—5s e 20s 


2 (s*x(s) sx(0) +(0))  sX(s) — x(0) +2X(s) = — - 


Using the initial values and solving for X(s), we find 


es os e 20s 


X(s) = —,——_... 
(s) (2s? + 5 +2) 

To determine the solution for x(f) we now need to find the inverse Laplace 
transform. The exponential functions can be dealt with using line 13 of Table 4.1. 
We write 

X(s) = (e-* —e-™) H(s), 


where i 
H(s) = ——~——_.. 
(s) s(2s* +s +2) 
Then using line 13, we have 
x(t) = us5(t)h(t — 5) — ug9(t)h(t — 20), (4.12) 


where h(t) = L~!{H(s)}. 
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To determine l(t) we need the partial fraction decomposition of H(s). Since the 
zeros of 2s* + s + 2 are complex, the correct decomposition is 


1 a bs+c 


s(2s*+s+2) 8 2s%+5+42’ 


yielding the equation 
1 = a(2s*+s4+2)+(bs+c)s; 


or after equating powers of s, 
2a+b=0, a+c=0, 2a=1, 


yielding a = 7 b= -—1,andc = —3. Therefore, 


172 s+5 
s 252 +542 


re! s+5 
— 2\s stttssi] 
2 
Inspecting Table 4.1, the first term can be transformed using line 2, and the second 


term can be transformed using lines 8 and 9, provided we complete the square of 
the denominator and then massage the numerator. That is, first we complete the 


square: 
ae re ee eae ae 
a aaa aie) 


H(s) = 


and next we write 


1 I 15 
s+5 (+2) + ay 
2 1 2 : 


Therefore, the function H(s) can be written as 


H(s) =3 1 (s+4) ( 1 ) iG 


s (s+z)?+R \vI5/ (84+49)2+2 


The first term is transformed using line 2, the second term using line 9 and the third 
term using line 8. We finally obtain 


h(t) = ; (1 ag (cos (V15t/4) + ag sin (V151/4)) ) ; (4.13) 


which when combined with (4.12) yields the rather complicated solution for x(t). 
We briefly comment that it is also possible to solve this example without using 
the Laplace transform. The key idea is that both x and x are continuous functions 
of t. Clearly from the form of the inhomogeneous term and the initial conditions, 
x(t) = 0 for 0 < t < 5. We then solve the ode between 5 < ft < 20 with the 
inhomogeneous term equal to unity and initial conditions x(5) = x(5) = 0. This 
requires first finding the general homogeneous solution, next finding a particular 
solution, and then adding the homogeneous and particular solutions and solving 
for the two unknown constants. To simplify the algebra, note that the best ansatz 
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to use to find the homogeneous solution is x(t) = e’'5), and not x(t) = e”. 


Finally, we solve the homogeneous ode for ¢ > 20 using as boundary conditions 
the previously determined values x(20) and x(20), where we have made use of the 
continuity of x and x. Here, the best ansatz to use is x(t) = e” (t-20) The student 
may benefit by trying this as an exercise and attempting to obtain a final solution 
that agrees with the form given by (4.12) and (4.13). 


Example: Solve 2% + x + 2x = 6(t —5) with x(0) = x(0) =0 


Here the inhomogeneous term is an impulse at time t = 5. Taking the Laplace 
transform of the ode using Table 4.1, and applying the initial conditions, 


2s* +5+2)X(s) =e, 
( 


so that 


The inverse Laplace transform may now be computed using lines 8 and 13 of Table 
4.1: 
2 : 
x(t) = agate (#-5)/4 sin (/15(t — 5) /4). (4.14) 

It is interesting to solve this example without using a Laplace transform. Clearly, 
x(t) = 0 up to the time of impulse at t = 5. Furthermore, after the impulse the ode 
is homogeneous and can be solved with standard methods. The only difficulty is 
determining the initial conditions of the homogeneous ode at t = 5t. 

When the inhomogeneous term is proportional to a delta-function, the solution 
x(t) is continuous across the delta-function singularity, but the derivative of the 
solution x(t) is discontinuous. If we integrate the second-order ode across the 
singularity at t = 5 and consider e — 0, the terms proportional to x and x on the 
left-hand side go to zero because of the continuity of x across the singularity, and 
only the second derivative term survives. Therefore, 


And as € + 0, we have x(5+) — x(5~) = 1/2. Since %(5~) = 0, the appropriate 
initial conditions immediately after the impulse force are x(5~) = 0 and x(57) = 
1/2. This result can be confirmed using (4.14). 
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Chapter 5 


Series solutions of 
homogeneous linear 
second-order differential 
equations 


Reference: Boyce and DiPrima, Chapter 5 


We consider the homogeneous linear second-order differential equation for y = 
y(x): 
P(x)y" + Q(x)y' + R(x)y = 0, (5.1) 


where P(x), Q(x) and R(x) are polynomials or convergent power series around 
x = Xo, with no common polynomial factors that could be divided out. The value 
xX = Xq is called an ordinary point of (5.1) if P(xo) € 0, and is called a singular point 
if P(x9) = 0. Singular points will later be further classified as regular singular points 
and irregular singular points. Our goal is to find two independent solutions of (5.1). 


5.1 Ordinary points 
If xo is an ordinary point of (5.1), then it is possible to determine two power series 


(ie., Taylor series) solutions for y = y(x) centered at x = xg. We illustrate the 
method of solution by solving two examples, with x9 = 0. 


Example: Find the general solution of y’ + y = 0. 
View tutorial on YouTube 
By now, you should know that the general solution is y(x) = ajcosx + a, sinx, 


with ag and a; constants. To find a power series solution about the point x9 = 0, we 
write 


y(x) = Do anx"; 
n=0 


and upon differentiating term-by-term 
co 
¥(x)= iz nayx" |, 
n=1 


and 


f(z) = y" n(n —1)ayx"-?. 


5.1. ORDINARY POINTS 


Substituting the power series for y and its derivatives into the differential equation 
to be solved, we obtain 


foe) 


Yo n(n —1)anx"* + Yo anx" = 0. (5.2) 
n=2 n=0 

The power-series solution method requires combining the two sums on the left- 

hand-side of (5.2) into a single power series in x. To shift the exponent of x2 in 

the first sum upward by two to obtain x”, we need to shift the summation index 

downward by two; that is, 


foe) lee} 
Yo n(n 1Vanx"? = Yo (n+2)(n + Lang2x". 
n=2 n=0 


We can then combine the two sums in (5.2) to obtain 


E ((n + 2)(n + 1)an424 an)" = 0. (5.3) 


n=0 


For (5.3) to be satisfied, the coefficient of each power of x must vanish separately. 
(This can be proved by setting x = 0 after successive differentiation.) We therefore 
obtain the recurrence relation 


m2 = 7 scl 2D ong: 


n+2)(n+1)’ 
We observe that even and odd coefficients decouple. We thus obtain two indepen- 
dent sequences starting with first term ag or a,. Developing these sequences, we 
have for the sequence beginning with ao: 


ao, 
a ay 
2= 2 0, 
a4 = ao = it a 
eo ee Tee 
ao : a4 = Das 
6 655 4 6! 0, 
and the general coefficient in this sequence for n = 0,1,2,... is 
—1)" 
a2n = Gyr 
Also, for the sequence beginning with a: 
M4, 
43=—- a 
3 3.0 
a5 = a3 = a 
a 2 ae De Ot, Ue ad 
1 1 
a7 = — 


7.60 = yi 
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and the general coefficient in this sequence for n = 0,1,2,... is 


a = ae a 


Using the principle of superposition, the general solution is therefore 


Ua 2 oa ca iy oe 
1 


n=0 


. i ae eS 
=a Be ae ee) ee | 


=agcosx + a, sinx, 


as expected. 
In our next example, we will solve the Airy’s Equation. This differential equation 
arises in the study of optics, fluid mechanics, and quantum mechanics. 


Example: Find the general solution of y” — xy = 0. 


View tutorial on YouTube 
With 
foe) 
= ‘Zz aax”, 
n=0 
the differential equation becomes 


foe} foe} 


Yo n(n —1)anx"? — Ye anx"** =0. (5.4) 


n=2 n=0 


We shift the first sum to x”+! by shifting the exponent up by three, i.e., 


CO co 
Yo n(n —1)anx"* = Y° (n+3)(n+2)an43x"*?. 
n=2 n=—-1 
When combining the two sums in (5.4), we separate out the extra n = —1 term in 


the first sum given by 2a. Therefore, (5.4) becomes 


2a2 + . ( + 3)(1+2)an+43 ay ja" =0. (5.5) 
n=0 


Setting coefficients of powers of x to zero, we first find az = 0, and then obtain the 


recursion relation ‘ 


An43 = (n+3)(n+2)" 


Three sequences of coefficients—those starting with either a9, a4; or a2,—decouple. 
In particular the three sequences are 


(5.6) 


A9,43,46,49,...-; 
Ay,44,A7,0409,---; 
42,45,4g,411..--- 
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Since a2 = 0, we find immediately for the last sequence 


ad as ag ay piss 0. 


We compute the first four nonzero terms in the power series with coefficients cor- 
responding to the first two sequences. Starting with aj, we have 


40, 

1 
43 3.200 
a6 = z a 
Gt §5o2 
ag = z ao; 
o Qepege he gea. 

and starting with a1, 

4, 

1 
Ha 
a7 = t a 
ss ee OT ead 

1 


M0 70005736 4a 


The general solution for y = y(x), can therefore be written as 


3x6 x9 gl. oh x10 
u(x) = a0 (1 “6 * 180 7 12960 * is) ra (: "12 7 504 7 45360 - =) 
= agyo(x) + a1y1(x). 


Suppose we would like to graph the solutions y = yo(x) and y = y;(x) versus x 
by solving the differential equation y” — xy = 0 numerically. What initial conditions 
should we use? Clearly, y = yo(x) solves the ode with initial values y(0) = 1 and 
y' (0) = 0, while y = y1(x) solves the ode with initial values y(0) = 0 and y'(0) = 1. 

The numerical solutions, obtained using MATLAB, are shown in Fig. 5.1. Note 
that the solutions oscillate for negative x and grow exponentially for positive x. 
This can be understood by recalling that y” + y = 0 has oscillatory sine and cosine 
solutions and y’’ — y = 0 has exponential hyperbolic sine and cosine solutions. 


5.2 Regular singular points: Cauchy-Euler equations 
View tutorial on YouTube 


The value x = x9 is called a regular singular point of the ode 


(x — x9)?y” + p(x)(x — xo)y/ + q(x)y = 0, (5.7) 


if p(x) and q(x) have convergent Taylor series about x = x9, ie., p(x) and q(x) can 
be written as a power-series in (x — xq): 


p(x) = pot pile — xp) + pa(x — x9) ++. 
q(x) = got qi(x — Xo) + qa(x— Xo) +--+, 
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Airy’s functions 


2 
-10 -8 -6 -4 -2 0 2 


Figure 5.1: Numerical solution of Airy’s equation. 


with py, and g, constants, and qg 4 0 so that (x — xg) is not a common factor 
of the coefficients. Any point x = x9 that is not an ordinary point or a regular 
singular point is called an irregular singular point. Many important differential 
equations of physical interest have regular singular points, and their solutions go 
by the generic name of special functions, with specific names associated with now 
famous mathematicians like Bessel, Legendre, Hermite, Laguerre and Chebyshev. 

Here, we will only consider the simplest ode with a regular singular point at 
x = 0. This ode is called a Cauchy-Euler equation, and has the form 


xy" + axy’ + By =0, (5.8) 


with a and B constants. Note that (5.7) reduces to a Cauchy-Euler equation (about 
xX = Xo) when one considers only the leading-order term in the Taylor series expan- 
sion of the functions p(x) and q(x). In fact, taking p(x) = po and q(x) = qo and 
solving the associated Cauchy-Euler equation results in at least one of the leading- 
order solutions to the more general ode (5.7). Often, this is sufficient to obtain initial 
conditions for numerical solution of the full ode. Students wishing to learn how to 
find the general solution of (5.7) can consult Boyce & DiPrima. 

An appropriate ansatz for (5.8) is y = x’, when x > 0 and y = (—x)" when 
x < 0, (or more generally, y = |x|" for all x), with r constant. After substitution into 
(5.8), we obtain for both positive and negative x 


r(r —1)|x|" + ar|x|" + Blx|" =0, 


and we observe that our ansatz is rewarded by cancelation of |x|". We thus obtain 
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the following quadratic equation for r: 
Pia rtp =0, (5.9) 


which can be solved using the quadratic formula. Three cases immediately appear: 
(i) real distinct roots, (ii) complex conjugate roots, (iii) repeated roots. Students may 
recall being in a similar situation when solving the second-order linear homoge- 
neous ode with constant coefficients. Indeed, it is possible to directly transform 
the Cauchy-Euler equation into an equation with constant coefficients so that our 
previous results can be used. 

The idea is to change variables so that the power law ansatz y = x” becomes an 
exponential ansatz. For x > 0, if we let x = e and y(x) = Y(é), then the ansatz 
y(x) = x" becomes the ansatz Y(é) = e’¢, appropriate if Y(€) satisfies a constant 
coefficient ode. If x < 0, then the appropriate transformation is x = —e$, since 
e& > 0. We need only consider x > 0 here and subsequently generalize our result 
by replacing x everywhere by its absolute value. 

We thus transform the differential equation (5.8) for y = y(x) into a differential 
equation for Y = Y(¢), using x = e, or equivalently, € = Inx. By the chain rule, 


dy _ ade 
dx dé dx 

_ lay 

~ x dé 
dY 

a2 BG 
e de’ 

so that symbolically, 

d _¢d 
—=e °_., 
dx dé 


The second derivative transforms as 


Py ed (AY 
dx2 dé dé 


pte fe ay 
° (e iz) 


Upon substitution of the derivatives of y into (5.8), and using x = e°, we obtain 


e% (e% (y" — ¥')) + ae (e*y’) + BY = Y" + (w—1)Y' + BY 
=0. 


As expected, the ode for Y = Y(@) has constant coefficients, and with Y = e’¢, the 
characteristic equation for the transformed differential equation is given by (5.9). We 
now directly transfer previous results obtained for the constant coefficient second- 
order linear homogeneous ode. 


5.2.1 Distinct real roots 


This simplest case needs no transformation. If (a — 1)? —4B > 0, then with r+ 
denoting the real roots of (5.9), the general solution is 


y(x) = ey|x|" + e|z|". 
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5.2.2 Distinct complex-conjugate roots 


If (a — 1)? — 4B < 0, we can write the complex roots of (5.9) as rz = A+ ip. Recall 
the general solution for Y = Y(¢) is given by 


Y(Z) =e (Acos é + Bsin pé) ; 
and upon transformation, and replacing x by |x|, 


y(x) = |x|*(Acos (In |x|) + Bsin (wIn|x|)). 


5.2.3 Repeated roots 


If (w — 1)? — 4B = 0, there is one real root r of (5.9). The general solution for Y is 


¥(Z) =e (cy +026), 
yielding 
y(x) = |x|" (cy + c2 In |x]). 


We now give examples illustrating these three cases. 


Example: Solve 2x*y + 3xy/ — y = 0 for 0 < x < 1 with two-point boundary 
condition y(0) = 0 and y(1) = 1. 
Since x > 0, we try y = x’ and obtain the characteristic equation 
0 =2r(r—1)+3r-1 
=2P°+r-1 
= (2r—1)(r+1). 
Since the characteristic equation has two real roots, the general solution is given by 


y(x) = cyx2 +egx 
We now encounter for the first time two-point boundary conditions, which can be 
used to determine the coefficients c,; and cy. Since y(0)=0, we must have cp = 0. 
Applying the remaining condition y(1) = 1, we obtain the unique solution 


Note that x = 0 is called a singular point of the ode since the general solution is 
singular at x = 0 when cp # 0. Our boundary condition imposes that y(x) is 
finite at x = 0 removing the singular solution. Nevertheless, y' remains singular at 
x = 0. Indeed, this is why we imposed a two-point boundary condition rather than 
specifying the value of y'(0) (which is infinite). 


Example: Solve x7y" + xy' + 7y = 0 with two-point boundary condition y(1) = 1 
and y(/e) =1. 
With the ansatz y = x", we obtain 

O=r(r—1)+r4+7° 


=P 4 72, 


CHAPTER 5. SERIES SOLUTIONS 73 


5.2. REGULAR SINGULAR POINTS: CAUCHY-EULER EQUATIONS 


so that r = ti. Therefore, with ¢ = Inx, we have Y(€) = Acos 7é + Bsin 7é, and 
the general solution for y(x) is 


y(x) = Acos(aInx) + Bsin(zInx). 


The first boundary condition y(1) = 1 yields A = 1. The second boundary condi- 
tion y(./e) = 1 yields B = 1. 


Example: Solve x*y" + 5xy’ + 4y = 0 with two-point boundary condition y(1) = 0 
and y(e) = 1. 
With the ansatz y = x", we obtain 
O=r(r—1)+5r+4 

=r+4r+4 

= ($2), 
so that there is a repeated root r = —2. With ¢ = Inx, we have Y(é) = (c1 + 
coé)e~*6, so that the general solution is 


cy tcp Inx 
x2 


y(x) = 


The first boundary condition y(1) = 0 yields c; = 0. The second boundary condi- 
tion y(e) = 1 yields cp = e*. The solution is therefore 


e2Inx 


y(x) = 


x2 
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Systems of coupled linear differential equations can result, for example, from lin- 
ear stability analyses of nonlinear equations, and from normal mode analyses of 
coupled oscillators. We will first consider the simplest case of a system of two 
coupled homogeneous linear first-order equations with constant coefficients. These 
two first-order equations are in fact equivalent to a single second-order equation, 
and the methods of Chapter 3 could be used for solution. Nevertheless, viewing 
the problem as a system of first-order equations introduces the important concept 
of the phase space, and can easily be generalized to higher-order linear systems. 
We will then discuss the physical problem of two coupled oscillators. 


6.1 Matrices, determinants and the eigenvalue prob- 
lem 


View a lecture on matrix addition and multiplication on YouTube 
View a lecture on determinants on YouTube 


We begin by reviewing some basic matrix algebra. A matrix with n rows and m 
columns is called an n-by-m matrix. Here, we need only consider the simple case 
of two-by-two matrices. 

A two-by-two matrix A, with two rows and two columns, can be written as 


a b 
Bie (° i) 
The first row has elements a and b, the second row has elements c and d. The first 
column has elements a and c; the second column has elements b and d. 


Matrices can be added and multiplied. Matrices can be added if they have the 
same dimension, and addition proceeds element by element, following 


a DV (ef \ fee Bae 

c ad gh) \ct+tg d+h)° 
Matrices can be multiplied if the number of columns of the left matrix equals the 
number of rows of the right matrix. A particular element in the resulting product 
matrix, say in row k and column I, is obtained by multiplying and summing the 


elements in row k of the left matrix with the elements in column / of the right 
matrix. For example, a two-by-two matrix can multiply a two-by-one column vector 


as follows 
a b\ (x\ _ fax+by 
c d}) \y)~ \cx+dy)° 
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The first row of the left matrix is multiplied against and summed with the first (and 
only) column of the right matrix to obtain the element in the first row and first 
column of the product matrix, and so on for the element in the second row and first 
column. The product of two two-by-two matrices is given by 


a b\ fe f\ _ faet+bg af+bh 
c d})\g h)” \ce+dg cf+dh)° 
A system of linear algebraic equations is easily represented in matrix form. For 


instance, with aj; and b; given numbers, and x; unknowns, the (two-by-two) system 
of equations given by 


Ay X1 + Ay2X2 = by, 


Ag X1 + A22X2 = by, 


can be written in matrix form as 
Ax =b, (6.1) 


A= & ean x= (es b= al (6.2) 
471 22 x2 bo 


When b = 0, we say that the system of equations given by (6.1) is homogeneous. 
We now ask the following question: When does there exist a nontrivial (not identi- 
cally zero) solution for x when the linear system is homogeneous? For the simplest 
case of a two-by-two matrix, we can try and solve directly the homogeneous linear 
system of equations given by 


where 


ax, + bxy = 0, 


cx; + dxy = 0. 6:3) 


Multiplying the first equation by d and the second by b, and subtracting the second 
equation from the first, results in 


(ad — bc)x, = 0. 


Similarly, multiplying the first equation by c and the second by a, and subtracting 
the first equation from the second, results in 


(ad — bc)xz = 0. 


Therefore, a nontrivial solution of (6.3) for x1 and x2 exists only if ad — be = 0. If 
(6.3) is expressed in matrix form and we define the determinant of the 2 x 2 matrix 


a b 
ro ‘ i) (6.4) 


to be det A = ad — bc, then we say that a nontrivial solution to (6.3) exists provided 
detA = 0. 
The same calculation may be repeated for a 3 x 3 matrix. If 


abe x41 
A = d e f 7 x= x2 , 
g hi x3 
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then there exists a nontrivial solution to Ax = 0 provided det A = 0, where det A = 
a(ei— fh) — b(di— fg) +c(dh—eg). The definition of the determinant can be further 
generalized to any n x n matrix, and is typically taught in a first course on linear 
algebra. 


We now consider the eigenvalue problem. For A an n X n matrix and v ann x 1 
column vector, the eigenvalue problem solves the equation 


Av = Av (6.5) 


for eigenvalues A; and corresponding eigenvectors v;. We rewrite the eigenvalue 
equation (6.5) as 


(A—Al)v =0, (6.6) 


where I is the n x n identity matrix, that is, the matrix with ones on the diagonal 
and zeros everywhere else. A nontrivial solution of (6.6) for A and v exists provided 


det (A — AI) =0. (6.7) 


Equation (6.7) is an n-th order polynomial equation in A, and is called the charac- 
teristic equation of A. The characteristic equation can be solved for the eigenvalues, 
and for each eigenvalue, a corresponding eigenvector can be determined directly 
from (6.5). 

We can demonstrate how to find the eigenvalues and eigenvectors of the 2 x 2 
matrix given by (6.4). We have 


0 = det (A — AI) 
a—A b | 


c d—xX 
= (a—A)(d— A) — be 
=A? —(a+d)A+ (ad — be). 


This characteristic equation can be more generally written as 
*—TrAA+detA =0, (6.8) 


where Tr A is the trace, or sum of the diagonal elements, of the matrix A. If A is an 
eigenvalue of A, then the corresponding eigenvector v may be found by solving 


a—xX b M1) 9 
c d—-A} \v.) ~~” 
where the equation of the second row will always be a multiple of the equation of 


the first row. The eigenvector v has arbitrary normalization, and we may always 
choose for convenience v; = 1. The equation from the first row is 


(a —A)v, + bv2 = 0, 


and with v; = 1, we find v2 = (A — a)/b. 
In the next section, we will see several examples of an eigenvector analysis. 
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6.2 Two coupled homogeneous linear first-order differ- 
ential equations 


We now consider the general system of differential equations given by 
Xy = ax, +bx, xX. = cx, + dx, (6.9) 
which can be written using matrix notation as 
x = Ax. (6.10) 


Before solving this system of odes using matrix techniques, I first want to show 
that we could actually solve these equations by converting the system into a single 
second-order equation. We take the derivative of the first equation and use both 
equations to write 


£1 = ax, + bry 

= ax, + b(cx, + dx2) 

= ax, + bex, + d(x, — ax) 
(a+ d)x, — (ad — bc)x}. 


The system of two first-order equations therefore becomes the following second- 
order equation: 
1 — (a+d)x1 + (ad — bc)x, = 0. 


If we had taken the derivative of the second equation instead, we would have ob- 
tained the identical equation for x2: 


¥y — (a+ d)xX_ + (ad — bc)x2 = 0. 


In general, a system of n first-order linear homogeneous equations can be converted 
into an equivalent n-th order linear homogeneous equation. Numerical methods 
usually require the conversion in reverse; that is, a conversion of an n-th order 
equation into a system of n first-order equations. 

With the ansatz x, = e*! or xp = e!, the second-order odes have the character- 
istic equation 


A? — (a+ d)A+ (ad — bc) =0. 


This is identical to the characteristic equation obtained for the matrix A from an 
eigenvalue analysis. 

We will see that (6.9) can in fact be solved by considering the eigenvalues and 
eigenvectors of the matrix A. We will demonstrate the solution for three separate 
cases: (i) eigenvalues of A are real and there are two linearly independent eigenvec- 
tors; (ii) eigenvalues of A are complex conjugates, and; (iii) A has only one linearly 
independent eigenvector. These three cases are analogous to the cases considered 
previously when solving the homogeneous linear constant-coefficient second-order 
equation. 


6.2.1 Distinct real eigenvalues 


We illustrate the solution method by example. 


78 CHAPTER 6. SYSTEMS OF EQUATIONS 


6.2. COUPLED FIRST-ORDER EQUATIONS 


Example: Find the general solution of %; = x1 + x2, X27 = 4x1 + x2. 


View tutorial on YouTube 


The equation to be solved may be rewritten in matrix form as 


Bix) fh 1) fa 
dt \ x2 ~A\4 1 xX)’ 
or using vector notation, written as (6.10). We take as our ansatz x(t) = ve, where 
v and A are independent of t. Upon substitution into (6.10), we obtain 
Ave# = Ave! ; 
and upon cancellation of the exponential, we obtain the eigenvalue problem 
Av = Av. (6.11) 


Finding the characteristic equation using (6.8), we have 


0 = det (A — AI) 


=A?-2A-3 
= (A—3)(A+1). 
Therefore, the two eigenvalues are Ay = 3 and Ay = —1. 


To determine the corresponding eigenvectors, we substitute the eigenvalues suc- 
cessively into 
(A—Al)v =0. (6.12) 


We will write the corresponding eigenvectors v1 and v2 using the matrix notation 


(v1 v2) _ ee 2) / 


021 «~(022 


where the components of vi and v2 are written with subscripts corresponding to 
the first and second columns of a 2 x 2 matrix. 

For A, = 3, and unknown eigenvector vy = (711 v21)!, where T denotes the 
transpose, we have from (6.12) 


—2011 + 021 = 0, 

4011 = 2021 = 0. 
Clearly, the second equation is just the first equation multiplied by —2, so only one 
equation is linearly independent. This will always be true, so for the 2 x 2 case 
we need only consider the first row of the matrix. The first eigenvector therefore 
satisfies v2; = 2011. Recall that an eigenvector is only unique up to multiplication 


by a constant: we may therefore take vj; = 1 for convenience. 
For Ay = —1, and unknown eigenvector v2 = (012 v22)', we have from (6.12) 


2012 + V2 = 0, 


so that V92 = —2v 2. Here, we take vj = 1. 
Therefore, our eigenvalues and eigenvectors are given by 


Ay = 3, n= (3); A,=—-1, v= (1). 
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Figure 6.1: Phase portrait with two real eigenvalues of opposite sign. 


Using the principle of superposition, the general solution to the ode is therefore 


x(t) = cyvy ce"! + Covner2t, 


or explicitly writing out the components, 


x(t) = cye* + ce, 


X2(t) = 2c,e% —2cne*. 


We can obtain a new perspective on the superposition of the two independent 
solutions by drawing a phase portrait, shown in Fig. 6.1, with “x-axis” x, and “y- 
axis” x9. Each curve corresponds to a different initial condition, and represents the 
trajectory of a hypothetical particle located at position (x1,x2) at different times t 
with velocity given by the differential equation. The dark lines represent trajectories 
along the direction of the eigenvectors. If cp = 0, the motion is along the eigenvector 
vi with x2 = 2x; and the motion with increasing time is away from the origin 
(arrows pointing out) since the eigenvalue A; = 3 > 0. If cj = 0, the motion is 
along the eigenvector v2 with x7 = —2x, and motion is towards the origin (arrows 
pointing in) since the eigenvalue A7 = —1 < 0. When the eigenvalues are real and 
of opposite signs, the origin is called a saddle point. Almost all trajectories (with 
the exception of those with initial conditions exactly satisfying x2(0) = —2x,(0)) 
eventually move away from the origin as t increases. 


Example: Find the general solution of x; = —3x, + V2x2, %2 = V2x1 — 2x2. 


The equations in matrix form are 


a(s)=(ya v2) (2): 
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Figure 6.2: Phase portrait with two real eigenvalues of same sign. 


The ansatz x = ve! leads to the eigenvalue problem Av = Av, with A the matrix 
above. The eigenvalues are determined from 
0 = det (A — AI) 
=AM+5A44 
= (A+4)(A+1). 


Therefore, the eigenvalues of A are A; = —4, Az = —1. Proceeding to determine 
the associated eigenvectors, for A; = —4, 


011 + V20y1 = 0; 


and for Az = —1, 
—2012 + V 2009 = 0. 


Taking the normalization vj; = 1 and vj2 = 1, we obtain for the eigenvalues and 
associated eigenvectors 


1 1 
m=-4=(_ 59): a= -1, v= (jp) 


The general solution to the ode is therefore 


(2) -o( sha) eeo( ae 


We show the phase portrait in Fig. 6.2. If cp = 0, the motion is along the 
eigenvector v1 with x. = —(/2/2)x, with eigenvalue A; = —4 < 0. If cy = 0, the 
motion is along the eigenvector v2 with x2 = J/2x1 with eigenvalue A2 = —1 < 0. 
When the eigenvalues are real and have the same sign, the origin is called a node. A 
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node may be attracting (as is the case here) or repelling depending on whether the 
eigenvalues are both negative or both positive. Observe that the trajectories collapse 
onto the v2 eigenvector since A; < Az < 0 and decay is more rapid along the v1 
direction. 


6.2.2 Distinct complex-conjugate eigenvalues 
Example: Find the general solution of x; = — 5X1 + x2, %2 = —xX1 — 5X2. 


View tutorial on YouTube 


The equations in matrix form are 


w(x) (4-4) @)- 


The ansatz x = ve"! leads to the equation 
0 = det (A — AI) 


5 
=M+A+-. 
+A Z 
Therefore, 1 = —1/2 +7; and we observe that the eigenvalues occur as a complex 
conjugate pair. We will denote the two eigenvalues as 


A=—5ti and A=—5-i. 


Now, for A a real matrix, if Av = Av, then AV = A¥. Therefore, the eigen- 
vectors also occur as a complex conjugate pair. The eigenvector v associated with 
eigenvalue A satisfies —iv, + v2 = 0, and normalizing with v; = 1, we have 


v=()), 


We have therefore determined two independent complex solutions to the ode, that 
is, 


ve and ve"! 


and we can form a linear combination of these two complex solutions to construct 
two independent real solutions. Namely, if the complex functions z(t) and Z(t) are 
written as 


z(t) = Re{z(t)} + iIm{z(t)}, Z(t) = Re{z(t)} — ilm{z(t)}, 
then two real functions can be constructed from the following linear combinations 
of z and z: 


Z 


2 +7 _Refz(t)} and <= Inte)}. 


2 
Thus the two real vector functions that can be constructed from our two complex 
vector functions are 


Re{ve*!} = Re { (;) ethan 
—l; 1 ‘oe 
=e 2 Re{ (;) (cost + isin’) } 
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Figure 6.3: Phase portrait with complex conjugate eigenvalues. 


and 


Im{ve*} = Him { (;) (cost + jsint) | 


Taking a linear superposition of these two real solutions yields the general solution 


to the ode, given by 
ae (4 ( ) ns: Gul 
—sint cos t 


The corresponding phase portrait is shown in Fig. 6.3. We say the origin is a 
spiral point. If the real part of the complex eigenvalue is negative, as is the case here, 
then solutions spiral into the origin. If the real part of the eigenvalue is positive, 
then solutions spiral out of the origin. 

The direction of the spiral—here, it is clockwise—can be determined easily. If 
we examine the ode with x; = 0 and x2 = 1, we see that x; = 1 and x) = —1/2. 
The trajectory at the point (0,1) is moving to the right and downward, and this 
is possible only if the spiral is clockwise. A counterclockwise trajectory would be 
moving to the left and downward. 


6.2.3 Repeated eigenvalues with one eigenvector 
Example: Find the general solution of xX; = x1 — x2, X2 = X1 + 3x2. 


View tutorial on YouTube 


The equations in matrix form are 


d (x 1 -1 x 
a(3)-G 3) (a): 9 
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The ansatz x = ve leads to the characteristic equation 
0 = det (A — AI) 
=N-4A+4 
= (A—2)?. 


Therefore, A = 2 is a repeated eigenvalue. The associated eigenvector is found from 
—v1 — 02 = 0, or v2 = —v,; and normalizing with v; = 1, we have 


We have thus found a single solution to the ode, given by 


m(t) =a (_1) 2% 


and we need to find the missing second solution to be able to satisfy the initial 
conditions. An ansatz of t times the first solution is tempting, but will fail. Here, we 
will cheat and find the missing second solution by solving the equivalent second- 
order, homogeneous, constant-coefficient differential equation. 

We already know that this second-order differential equation for x;(t) has a 
characteristic equation with a degenerate eigenvalue given by A = 2. Therefore, the 
general solution for x; is given by 


x1(t) = (cy + teg)e”*. 
Since from the first differential equation, x2 = x, — x1, we compute 
1 = (2cy + (1+ 2#)c2)e*, 
so that 


xg =X - x4 
= (cy + tez)e* — (2c, + (1+ 2t)c2)e”* 


= —ee** + (—1— te”. 


Combining our results for x; and x2, we have therefore found 


(8) <a(A)e+e[(9)+(A)e 


Our missing linearly independent solution is thus determined to be 


x(t) = co (e& + & ] ae (6.14) 


The second term of (6.14) is just ¢ times the first solution; however, this is not 
sufficient. Indeed, the correct ansatz to find the second solution directly is given by 


x = (w+ tv) aie (6.15) 
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Figure 6.4: Phase portrait with only one eigenvector. 


where A and v are the eigenvalue and eigenvector of the first solution, and w is an 
unknown vector to be determined. To illustrate this direct method, we substitute 
(6.15) into x = Ax, assuming Av = Av . Canceling the exponential, we obtain 


v+A(w+tv) = Aw4 Atv. 
Further canceling the common term Atv and rewriting yields 
(A—AI)w=v. (6.16) 


If A has only a single linearly independent eigenvector v, then (6.16) can be solved 
for w (otherwise, it cannot). Using A, A and v of our present example, (6.16) is the 
system of equations given by 


—1 —-1)\ /w,\ _ 1 
1 1) \wo)  \-1)° 
The first and second equation are the same, so that w2 = —(w, + 1). Therefore, 
w= A 
“=i 


=e()+(4), 


Notice that the first term repeats the first found solution, i-e., a constant times the 
eigenvector, and the second term is new. We therefore take w, = 0 and obtain 


— 0 
=( 4] 7 
as before. 


The phase portrait for this ode is shown in Fig. 6.4. The dark line is the single 
eigenvector v of the matrix A. When there is only a single eigenvector, the origin is 
called an improper node. 
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6.3. Normal modes 
View tutorials on YouTube: Part 1 Part 2 


We now consider an application of the eigenvector analysis to the coupled mass- 
spring system shown in Fig. 6.5. The position variables x; and x2 are measured 
from the equilibrium positions of the masses. Hooke’s law states that the spring 
force is linearly proportional to the extension length of the spring, measured from 
equilibrium. By considering the extension of the spring and the sign of the force, 
we write Newton’s law F = ma separately for each mass: 


mx = —kx1 = K(x4 = x2), 
mxXq = —kx — K(x2 — x1). 


Further rewriting by collecting terms proportional to x; and x2 yields 


mx, = —(k+ K)x, + Kxo, 
mXp = Kx, = (k+ K)x». 


The equations for the coupled mass-spring system form a system of two second- 
order linear homogeneous odes. In matrix form, mk = Ax, or explicitly, 


WB (2)-(82" ake) (0 


In analogy to a system of first-order equations, we try the ansatz x = ve! and upon 


substitution into (6.17) we obtain the eigenvalue problem Av = Av, with A = mr’. 


The eigenvalues are determined by solving the characteristic equation 
0 = det (A — AI) 
_ |-(k+K)—-A K 
a K —(k+K)—-A 
= (A+k+K) — R?. 


The solution for A is 
A=-k-KxK, 


and the two eigenvalues are 


Ay =—k, Ay =—(k+2K). 


The corresponding values of r in our ansatz x = ve" f with r = +V/A/m, are 


ry =ivk/m, 7, 2 =i\/(kK+2K)/m, 7p. 


Since the values of r are pure imaginary, we know that x(t) and x(t) will oscillate 
with angular frequencies w, = Im{r,} and w2 = Im{rp}, that is, 


wy =Vk/m, wo = 1/(k+2K)/m. 


The positions of the oscillating masses in general contain time dependencies of the 
form sin wt, cos wy t, and sin wt, cos Wot. 
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Figure 6.5: Coupled harmonic oscillators. 


It is of further interest to determine the eigenvectors, or so-called normal modes 
of oscillation, associated with the two distinct angular frequencies. With specific 
initial conditions proportional to an eigenvector, the mass will oscillate with a single 
frequency. The eigenvector with eigenvalue A, satisfies 


—Koy + Koy, = 0, 


so that v1; = v2. The normal mode with frequency w,; = Vk/m thus follows a 
motion where x; = x2. Referring to Fig. 6.5, during this motion the center spring 
length does not change, which is why the frequency of oscillation is independent 
of K. 

Next, we determine the eigenvector with eigenvalue A: 


Koz, + Koz) = 0, 


so that v2] = —v22. The normal mode with frequency w2 = \/(k+2K)/m thus 


follows a motion where x; = —x2. Again referring to Fig. 6.5, during this motion 
the two equal masses symmetrically push or pull against each side of the middle 
spring. 


A general solution for x(t) can be constructed from the eigenvalues and eigen- 
vectors. Our ansatz was x = ve", and for each of two eigenvectors v, we have a 
pair of complex conjugate values for r. Accordingly, we first apply the principle of 
superposition to obtain four real solutions, and then apply the principle again to 
obtain the general solution. With w, = Vk/m and wy = \/(k+2K)/m, the general 
solution is given by 


(3) = (;) (A cos w t+ Bsinwyt) + &) (Ccos wat + Dsinwyt), 


where the now real constants A, B, C, and D can be determined from the four 
independent initial conditions, x1(0), x2(0), %1(0), and x2(0). 
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Nonlinear differential 
equations and bifurcation 
theory 


Reference: Strogatz, Sections 2.2, 2.4, 3.1, 3.2, 3.4, 6.3, 6.4, 8.2 


We now turn our attention to nonlinear differential equations. In particular, we 
study how small changes in the parameters of a system can result in qualitative 
changes in the dynamics. These qualitative changes in the dynamics are called 
bifurcations. To understand bifurcations, we first need to understand the concepts 
of fixed points and stability. 


7.1 Fixed points and stability 


7.1.1 One dimension 


View tutorial on YouTube 


Recall the non-linear logistic equation, x = x(x — 1), considered in Section 2.4.6. 
This is an example of the general one-dimensional differential equation for x = x(f) 
given by 

x = f(x). (7.1) 
We say that x. is a fixed point, or equilibrium point, of (7.1) if f(x.) = 0. At the fixed 
point, x = 0. The terminology fixed point is used since the solution to (7.1) with 
initial condition x(0) = x, is x(t) = x, for all time t. 

A fixed point, however, can be stable or unstable. A fixed point is said to be 
stable if a small perturbation of the solution from the fixed point decays in time; it is 
said to be unstable if a small perturbation grows in time. We can determine stability 
by a linear analysis. Let x = x» + e(t), where € represents a small perturbation of 
the solution from the fixed point x,. Because x, is a constant, x = €; and because x, 
is a fixed point, f(x,) = 0. Taylor series expanding the function f about the fixed 
point x = x, results in the differential equation 


é= f(x. +e) 
= f(x*) tefl (ts) +... 
SEF ie ie 
The omitted terms in the Taylor series are proportional to e?, and can be made 
negligible over a short time interval with respect to the kept term, proportional to 
€, by taking the initial perturbation, e(0), sufficiently small. Therefore, at least over 


short times, the differential equation to be considered, é = f’(x.)e, is linear and has 
by now the familiar solution 


e(t) = €(O)ef Ow, 
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The exponential function decays exponentially if f(x.) < 0, and we say the fixed 
point is stable. The exponential function grows exponentially if f' "(x,) > 0, and we 
say the fixed point is unstable. If f’(x..) = 0, we say the fixed point is marginally 
stable, and the next higher-order term in the Taylor series must be considered. 


Example: Find all the fixed points of the logistic equation x = x(1— x) and 
determine their stability. 


There are two fixed points at which x = 0, given by x, = 0 and x, = 1. Stability 
of these equilibrium points may be determined by considering the derivative of 
f(x) = x(1 — x). We have f(x) = 1 - 2x. Therefore, f’(0) = 1 > 0 so that x, = 0 is 
an unstable fixed point, and f’(1) = —1 < 0 so that x, = 1 is a stable fixed point. 
Indeed, we have previously found that all solutions approach the stable fixed point 
asymptotically. 


7.1.2 Two dimensions 


View tutorial on YouTube 


The idea of fixed points and stability can be extended to higher-order systems of 
odes. Here, we consider a two-dimensional system and will need to make use of 
the two-dimensional Taylor series of a function F(x,y). In general, the Taylor series 
of F(x,y) expanded about the origin is given by 


OF OF 5 (#55 o°F et) 


ox Yay FON” oa2 r Yay ay TY oy 


| 
Te ost asety: 


F(x,y) =F+x 


where the function F and all of its partial derivatives on the right-hand side are 
evaluated at the origin. Note that the Taylor series is constructed so that all the 
partial derivatives of the left-hand side evaluated at the origin match those of the 
right-hand side. 

We now consider the two-dimensional system given by 


x= f(x,y), Y= (x,y). (7.2) 


The point (x, ¥,) is said to be a fixed point of (7.2) if f(x», ys.) = Oand 9(xx, y+) = 0. 
Again, the local stability of a fixed point can be determined by a linear analysis. 
We let x(t) = x. +e(t) and y(t) = y, + 6(t), where € and 6 are small independent 
perturbations from the fixed point. Making use of the two dimensional Taylor series 
of f(x,y) and g(x,y) about the fixed point, or equivalently about (€,6) = (0,0), we 
have 


€ = f (Xs +€,Yx +0) 


See 
ftes ay bee. 


x 


Pee i 


b = o(x« + €,Ys +4) 
Og 08 
Poy, te 


=e Fea, 


= 208 5 498 
=a Se gg Per 
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Figure 7.1: Phase portrait for two-dimensional nonlinear system. 


where in the Taylor series f, g and all their partial derivatives are evaluated at the 
fixed point (x.,y.). Neglecting higher-order terms in the Taylor series, we thus 
have a system of odes for the perturbation, given in matrix form as 


d fe oo FN re 
a(5)-(% #) (5): ve) 
ox 


The two-by-two matrix in (7.3) is called the Jacobian matrix at the fixed point. An 
eigenvalue analysis of the Jacobian matrix will typically yield two eigenvalues A; 
and Az. These eigenvalues may be real and distinct, complex conjugate pairs, or 
repeated. The fixed point is stable (all perturbations decay exponentially) if both 
eigenvalues have negative real parts. The fixed point is unstable (some perturba- 
tions grow exponentially) if at least one of the eigenvalues has a positive real part. 
Fixed points can be further classified as stable or unstable nodes, unstable saddle 
points, stable or unstable spiral points, or stable or unstable improper nodes. 


3 


Example: Find all the fixed points of the nonlinear system x = x(3 — x — 2y), 
y = y(2—x-—y), and determine their stability. 


View tutorial on YouTube 


The fixed points are determined by solving 


f(xy) =x3—x—-2y)=0, guy) =y2—x—y) =0. 
Evidently, (x,y) = (0,0) is a fixed point. On the one hand, if only x = 0, then 
the equation g(x,y) = 0 yields y = 2. On the other hand, if only y = 0, then the 
equation f(x,y) = 0 yields x = 3. If both x and y are nonzero, then we must solve 
the linear system 
x+2y=3, x+y=2, 
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x (1, 
the four fixed points (x.,y.) = (0,0), (0,2), (3,0), (1,1). The Jacobian matrix is 


given by 
of of 
dx oy \ _ ( 3—2x—2y —2x 
og og J -y 2-x-2y }* 


Stability of the fixed points may be considered in turn. With J, the Jacobian matrix 
evaluated at the fixed point, we have 


(0) = 00): K=(5 9): 


The eigenvalues of J, are A = 3,2 so that the fixed point (0,0) is an unstable node. 
Next, 


(#46) = 0.2): =( 2 _3 )- 


The eigenvalues of J. are A = —1,—2 so that the fixed point (0,2) is a stable node. 
Next, 


(4.6) =@.0): = (GIy )- 


The eigenvalues of J,, are A = —3,—1 so that the fixed point (3,0) is also a stable 
node. Finally, 


(0) =): =( 2) Tf). 


The characteristic equation of J, is given by (—1— A)* —2 = 0, so that A = -14+ V2. 
Since one eigenvalue is negative and the other positive the fixed point (1,1) is an 
unstable saddle point. From our analysis of the fixed points, one can expect that all 
solutions will asymptote to one of the stable fixed points (0,2) or (3,0), depending 
on the initial conditions. 

It is of interest to sketch the phase portrait for this nonlinear system. The eigen- 
vectors associated with the unstable saddle point (1,1) determine the directions of 
the flow into and away from this fixed point. The eigenvector associated with the 
positive eigenvalue A; = —1+ /2 can be determined from the first equation of 
(Jx — AqI)v1 = 0, or 


~V 20, — 2012 = 0, 


so that vy. = — (/2 /2)v11. The eigenvector associated with the negative eigenvalue 
Ay = —-1— v2 satisfies 022 = (/2/2)v21. The eigenvectors give the slope of the 
lines with origin at the fixed point for incoming (negative eigenvalue) and outgo- 
ing (positive eigenvalue) trajectories. The outgoing trajectories have negative slope 
—1/2/2 and the incoming trajectories have positive slope J2/2. A rough sketch of 
the phase portrait can be made by hand (as demonstrated in class). Here, a com- 
puter generated plot obtained from numerical solution of the nonlinear coupled 
odes is presented in Fig. 7.1. The curve starting from the origin and at infinity, 
and terminating at the unstable saddle point is called the separatrix. This curve 
separates the phase space into two regions: initial conditions for which the solution 
asymptotes to the fixed point (0,2), and initial conditions for which the solution 
asymptotes to the fixed point (3,0). 
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(a) dx/dt 
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r<0 r=0 r>0 
(:) ied te = 


Figure 7.2: Saddlenode bifurcation. (a) x versus x; (b) bifurcation diagram. 


7.2 One-dimensional bifurcations 


A bifurcation occurs in a nonlinear differential equation when a small change in 
a parameter results in a qualitative change in the long-time solution. Examples of 
bifurcations are when fixed points are created or destroyed, or change their stability. 

We now consider four classic bifurcations of one-dimensional nonlinear differen- 
tial equations: saddle-node bifurcation, transcritical bifurcation, supercritical pitch- 
fork bifurcation, and subcritical pitchfork bifurcation. The corresponding differen- 
tial equation will be written as 

& = fr(x), 


where the subscript r represents a parameter that results in a bifurcation when 
varied across zero. The simplest differential equations that exhibit these bifurcations 
are called the normal forms, and correspond to a local analysis (i.e., Taylor series) of 
more general differential equations around the fixed point, together with a possible 
rescaling of x. 


7.2.1 Saddle-node bifurcation 


View tutorial on YouTube 


The saddle-node bifurcation results in fixed points being created or destroyed. The 
normal form for a saddle-node bifurcation is given by 


x=rt+x. 


The fixed points are x, = +,/—r. Clearly, two real fixed points exist when r < 0 
and no real fixed points exist when r > 0. The stability of the fixed points when 
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r<0 r=0 r>0 


Figure 7.3: Transcritical bifurcation. (a) x versus x; (b) bifurcation diagram. 


r < 0 are determined by the derivative of f(x) = r+x?, given by f’(x) = 2x. 
Therefore, the negative fixed point is stable and the positive fixed point is unstable. 

Graphically, we can illustrate this bifurcation in two ways. First, in Fig. 7.2(a), 
we plot x versus x for the three parameter values corresponding to r < 0, r = 0 and 
r > 0. The values at which x = 0 correspond to the fixed points, and arrows are 
drawn indicating how the solution x(t) evolves (to the right if x > 0 and to the left 
if x < 0). The stable fixed point is indicated by a filled circle and the unstable fixed 
point by an open circle. Note that when r = 0, solutions converge to the origin from 
the left, but diverge from the origin on the right. Second, in Fig. 7.2(b), we plot a 
bifurcation diagram illustrating the fixed point x, versus the bifurcation parameter 
r. The stable fixed point is denoted by a solid line and the unstable fixed point by 
a dashed line. Note that the two fixed points collide and annihilate at r = 0, and 
there are no fixed points for r > 0. 


7.2.2 Transcritical bifurcation 


View tutorial on YouTube 


A transcritical bifurcation occurs when there is an exchange of stabilities between 
two fixed points. The normal form for a transcritical bifurcation is given by 


x= rx — x. 


The fixed points are x, = 0 and x, = r. The derivative of the right-hand-side is 


f'(x) = r— 2x, so that f’(0) = r and f'(r) = —r. Therefore, for r < 0, x. = 0 
is stable and x, = r is unstable, while for r > 0, x, = r is stable and x, = 0 is 
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Figure 7.4: Supercritical pitchfork bifurcation. (a) x versus x; (b) bifurcation diagram. 


unstable. The two fixed points thus exchange stability as r passes through zero. 
The transcritical bifurcation is illustrated in Fig. 7.3. 


7.2.3 Supercritical pitchfork bifurcation 


View tutorial on YouTube 


The pitchfork bifurcations occur in physical models where fixed points appear and 
disappear in pairs due to some intrinsic symmetry of the problem. Pitchfork bi- 
furcations can come in one of two types. In the supercritical bifurcation, a pair of 
stable fixed points are created at the bifurcation (or critical) point and exist after 
(super) the bifurcation. In the subcritical bifurcation, a pair of unstable fixed points 
are created at the bifurcation point and exist before (sub) the bifurcation. 

The normal form for the supercritical pitchfork bifurcation is given by 


x=rx—x°. 


Note that the linear term results in exponential growth when r > 0 and the non- 
linear term stabilizes this growth. The fixed points are x, = 0 and x, = + Jr, the 
latter fixed points existing only when r > 0. The derivative of f is f’(x) = r — 3x? 
so that f’(0) = rand f’(+./r) = —2r. Therefore, the fixed point x. = 0 is stable for 
r <0 and unstable for r > 0 while the fixed points x = +,/r exist and are stable for 
r > 0. Notice that the fixed point x, = 0 becomes unstable as r crosses zero and two 
new stable fixed points x, = +,/r are born. The supercritical pitchfork bifurcation 
is illustrated in Fig. 7.4. 
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xX, 


Figure 7.5: Subcritical pitchfork bifurcation. 


7.2.4 Subcritical pitchfork bifurcation 
View tutorial on YouTube 
In the subcritical case, the cubic term is destabilizing. The normal form (to order 
x3) is 
k= rx+x°. 


The fixed points are x, = 0 and x, = +,/—,, the latter fixed points existing only 
when r < 0. The derivative of the right-hand-side is f’(x) = r+ 3x? so that f’(0) =r 
and f’(+./—r) = —2r. Therefore, the fixed point x, = 0 is stable for r < 0 and 
unstable for r > 0 while the fixed points x = +./—r exist and are unstable for 
r <0. There are no stable fixed points when r > 0. 

The absence of stable fixed points for r > 0 indicates that the neglect of terms of 
higher-order in x than x3 in the normal form may be unwarranted. Keeping to the 
intrinsic symmetry of the equations (only odd powers of x) we can add a stabilizing 
nonlinear term proportional to x°. The extended normal form (to order x°) is 


Y=rx+x2—x, 


and is somewhat more difficult to analyze. The fixed points are solutions of 
meg a at, 


The fixed point x, = 0 exists for all r, and four additional fixed points can be found 


from the solutions of the quadratic equation in x: 
1 
te = 44/501 VIT®). 


These fixed points exist only if x, is real. Clearly, for the inner square-root to be 
real, r > —1/4. Also observe that 1 — 1+ 4r becomes negative for r > 0. We thus 
have three intervals in r to consider, and these regions and their fixed points are 


r< -j : X, =0 (one fixed point); 
1 1 er: : 
gl 0: x»%=0, %=+ 5 (1 +V1+4r) (five fixed points); 
r>0) x.=0, x= x(t 1+4r) (three fixed points). 
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Stability is determined from f’(x) = r+ 3x* —5x*. Now, f’(0) = r so x. = 0 is 
stable for r < 0 and unstable for r > 0. The calculation for the other four roots can 
be simplified by noting that x, satisfies r+ x2 — x4 = 0, or xt =r + x2. Therefore, 
f(4)=artat =55 
=7+3x2 —50 +22) 
= —4y — 2x2 
S20), 


With x2 = 5(1+ 1+ -4r), we have 


Puy (2r+ atid viv) 
ms VI+%) 


2 
4r) + 
VI+4 ( VI+4r+1). 


Clearly, the plus root is always stable since f(x.) < 0. The minus root exists only 
for —j <r < and is unstable since f’(x..) > 0. We summarize the stability of the 
various fixed points: 


= 


r< -; : XxX, =0 (stable); 


-; <r<0: x,=0, (stable) 
Xia (1 +/1+4r) (stable); 
x, = (1 v1+4r) (unstable); 


r>0O: x. =0 (unstable) 


1 
Xe = 5 +V1+4r) (stable). 


The bifurcation diagram is shown in Fig. 7.5, and consists of a subcritical pitch- 
fork bifurcation at r = 0 and two saddle-node bifurcations at r = —1/4. We can 
imagine what happens to x as r increases from negative values, supposing there is 
some small noise in the system so that x = x(t) will diverge from unstable fixed 
points. For r < —1/4, the equilibrium value of x is x, = 0. As r increases into 
the range —1/4 < r < 0, x will remain at x, = 0. However, a catastrophe occurs 
as soon as r > 0. The x, = 0 fixed point becomes unstable and the solution will 
jump up (or down) to the only remaining stable fixed point. Such behavior is called 
a jump bifurcation. A similar catastrophe can happen as r decreases from positive 
values. In this case, the jump occurs as soon as r < —1/4. 

Since the stable equilibrium value of x depends on whether we are increasing or 
decreasing 1, we say that the system exhibits hysteresis. The existence of a subcriti- 
cal pitchfork bifurcation can be very dangerous in engineering applications since a 
small change in a problem’s parameters can result in a large change in the equilib- 
rium state. Physically, this can correspond to a collapse of a structure, or the failure 
of a component. 
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7.2.5 Application: a mathematical model of a fishery 


View tutorial on YouTube 


We illustrate the utility of bifurcation theory by analyzing a simple model of a fish- 
ery. We utilize the logistic equation (see §2.4.6) to model a fish population in the 
absence of fishing. To model fishing, we assume that the government has estab- 
lished fishing quotas so that at most a total of C fish per year may be caught. We 
assume that when the fish population is at the carrying capacity of the environ- 
ment, fisherman can catch nearly their full quota. When the fish population drops 
to lower values, fish may be harder to find and the catch rate may fall below C, 
eventually going to zero as the fish population diminishes. Combining the logistic 
equation together with a simple model of fishing, we propose the mathematical 
model 


(7.4) 


where N is the fish population size, t is time, r is the maximum potential growth 
rate of the fish population, K is the carrying capacity of the environment, C is the 
maximum rate at which fish can be caught, and A is a constant satisfying A < K 
that is used to model the idea that fish become harder to catch when scarce. 

We nondimensionalize (7.4) using x = N/K, tT = rt,c = C/rK, « = A/K: 


(7.5) 


Note that 0 < x < 1,c > 0 and 0 < a < 1. We wish to qualitatively describe the 
equilibrium solutions of (7.5) and the bifurcations that may occur as the nondimen- 
sional catch rate c increases from zero. Practically, a government would like to issue 
each year as large a catch quota as possible without adversely affecting the number 
of fish that may be caught in subsequent years. 

The fixed points of (7.5) are x, = 0, valid for all c, and the solutions to x? — (1 — 


a)x+(c—«) =0,or 
[aa a) 25 4/ (1+ a)? —4c| . (7.6) 


The fixed points given by (7.6) are real only if c < i +«)*. Furthermore, the 
minus root is greater than zero only if c > a. We therefore need to consider three 
intervals over which the following fixed points exist: 


XxX, = 


Nir 


O<c<a: x,=0, r= 500 a) (1+«)2 ae]; 


pee ae SO. m= 5 /G-a)+ (1+ 0)? —4e] 


x = 0) 


The stability of the fixed points can be determined with rigor analytically or graph- 
ically. Here, we simply apply biological intuition together with knowledge of the 
types of one dimensional bifurcations. An intuitive argument is made simpler if we 
consider c decreasing from large values. When c is large, that is c > x(1 + a)*, too 
many fish are being caught and our intuition suggests that the fish population goes 
extinct. Therefore, in this interval, the single fixed point x, = 0 must be stable. As 
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Figure 7.6: Fishery bifurcation diagram. 


c decreases, a bifurcation occurs at c = 7(1 + «)? introducing two additional fixed 
points at x, = (1—.«)/2. The type of one dimensional bifurcation in which two 
fixed points are created as a square root becomes real is a saddlenode bifurcation, 
and one of the fixed points will be stable and the other unstable. Following these 
fixed points as c —> 0, we observe that the plus root goes to one, which is the appro- 
priate stable fixed point when there is no fishing. We therefore conclude that the 
plus root is stable and the minus root is unstable. As c decreases further from this 
bifurcation, the minus root collides with the fixed point x, = 0 at c = a. This ap- 
pears to be a transcritical bifurcation and assuming an exchange of stability occurs, 
we must have the fixed point x, = 0 becoming unstable for c < a. The resulting 
bifurcation diagram is shown in Fig. 7.6. 

The purpose of simple mathematical models applied to complex ecological prob- 
lems is to offer some insight. Here, we have learned that overfishing (in the model 
c> x(1 + )*) during one year can potentially result in a sudden collapse of the 
fish catch in subsequent years, so that governments need to be particularly cautious 
when contemplating increases in fishing quotas. 


7.3 Two-dimensional bifurcations 


All the one-dimensional bifurcations can also occur in two-dimensions along one 
of the directions. In addition, a new type of bifurcation can also occur in two- 
dimensions. Suppose there is some control parameter p. Furthermore, suppose 
that for 1 < 0, a two-dimensional system approaches a fixed point by exponentially- 
damped oscillations. We know that the Jacobian matrix at the fixed point with p < 0 
will have complex conjugate eigenvalues with negative real parts. Now suppose 
that when pt > 0 the real parts of the eigenvalues become positive so that the 
fixed point becomes unstable. This change in stability of the fixed point is called 
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a Hopf bifurcation. The Hopf bifurcation comes in two types: supercritical Hopf 
bifurcation and subcritical Hopf bifurcation. For the supercritical Hopf bifurcation, 
as pt increases slightly above zero, the resulting oscillation around the now unstable 
fixed point is quickly stabilized at small amplitude. This stable orbit is called a limit 
cycle. For the subcritical Hopf bifurcation, as p increases slightly above zero, the 
limit cycle immediately jumps to large amplitude. 


7.3.1 Supercritical Hopf bifurcation 


A simple example of a supercritical Hopf bifurcation can be given in polar coordi- 
nates: 


ram, 


6=w-+bdr, 


where x = rcos@ and y = rsin@. The parameter pt controls the stability of the fixed 
point at the origin, the parameter w is the frequency of oscillation near the origin, 
and the parameter b determines the dependence of the oscillation frequency at 
larger amplitude oscillations. Although we include b for generality, our qualitative 
analysis of these equations will be independent of b. 


The equation for the radius is of the form of the supercritical pitchfork bifurca- 
tion. The fixed points are r, = 0 and r, = JE (note that r > 0), and the former 
fixed point is stable for p < 0 and the latter is stable for pp > 0. The transition 
of the eigenvalues of the Jacobian from negative real part to positive real part can 
be seen if we transform these equations to cartesian coordinates. We have using 
pe a gh. 2, 

x = ¢cos6 — 6rsiné 
= (ur—r°) cosé — (w + br*)rsin®é 
= px (x + y?)x— wy bx? +y)y 
= px — wy — (x7 +y?)(x + by); 

y =7sin0 + 6rcosé 
= (ur —r?) sin@ + (w + br*)rcosé 
= py (27? + )y + wx + W(x? +y°)x 
= wx + py — (2° +y°)(y — bx). 


The stability of the origin is determined by the Jacobian matrix evaluated at the 
origin. The nonlinear terms in the equation vanish and the Jacobian matrix at the 


origin is given by 
= how 
a ( woo ). 


The eigenvalues are the solutions of (uw — A)? + w? = 0, or A = ptiw. As p in- 
creases from negative to positive values, exponentially damped oscillations change 
into exponentially growing oscillations. The nonlinear terms in the equations stabi- 
lize the growing oscillations into a limit cycle. 
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7.3.2 Subcritical Hopf bifurcation 
The analogous example of a subcritical Hopf bifurcation is given by 


f= ={P, 


6=w + br’. 


Here, the equation for the radius is of the form of the subcritical pitchfork bifurca- 
tion. As p increases from negative to positive values, the origin becomes unstable 
and exponentially growing oscillations increase until the radius reaches a stable 
fixed point far away from the origin. In practice, it may be difficult to tell ana- 
lytically if a Hopf bifurcation is supercritical or subcritical from the equations of 
motion. Computational solution, however, can quickly distinguish between the two 


types. 
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Chapter 8 


Partial differential equations 


Reference: Boyce and DiPrima, Chapter 10 


Differential equations containing partial derivatives with two or more independent 
variables are called partial differential equations (pdes). These equations are of 
fundamental scientific interest but are substantially more difficult to solve, both 
analytically and computationally, than odes. 

In this chapter, we begin by deriving two fundamental pdes: the diffusion equa- 
tion and the wave equation, and show how to solve them with prescribed boundary 
conditions using the technique of separation of variables. We then discuss solutions 
of the two-dimensional Laplace equation in cartesian and polar coordinates, and 
finish with a lengthy discussion of the Schrédinger equation, a partial differential 
equation fundamental to both physics and chemistry. 


8.1 Derivation of the diffusion equation 


To derive the diffusion equation in one spacial dimension, we imagine a still liquid 
in a long pipe of constant cross sectional area. A small quantity of dye is placed in 
a cross section of the pipe and allowed to diffuse up and down the pipe. The dye 
diffuses from regions of higher concentration to regions of lower concentration. 

We define u(x,t) to be the concentration of the dye at position x along the 
pipe, and we wish to find the pde satisfied by u. It is useful to keep track of 
the units of the various quantities involved in the derivation and we introduce the 
bracket notation [X] to mean the units of X. Relevant dimensional units used in the 
derivation of the diffusion equation are mass m, length I, and time t. Assuming that 
the dye concentration is uniform in every cross section of the pipe, the dimensions 
of concentration used here are [u] = m/l. 

The mass of dye in the infinitesimal pipe volume located between position x1 
and position x2 at time t, with x1 < x < x2, is given to order Ax = x2 — x1 by 


M = u(x,t)Ax. 


The mass of dye in this infinitesimal pipe volume changes by diffusion into and out 
of the cross sectional ends situated at position x; and x2 (Figure 8.1). We assume 
the rate of diffusion is proportional to the concentration gradient, a relationship 
known as Fick’s law of diffusion. Fick’s law of diffusion assumes the mass flux J, 
with units [J] = m/t across a cross section of the pipe is given by 


J=—Dux, (8.1) 
where the diffusion constant D > 0 has units [D] = I*/t, and we have used the 
notation ux = du/dx. The mass flux is opposite in sign to the gradient of concen- 
tration so that the flux is from high concentration to low concentration. The time 
rate of change in the mass of dye between x and x2 is given by the difference be- 
tween the mass flux into and the mass flux out of the infinitesimal cross sectional 


103 


8.2. DERIVATION OF THE WAVE EQUATION 


xX X 


Figure 8.1: Derivation of the diffusion equation. 


volume. A positive mass flux signifies diffusion from left to right. Therefore, the 
time rate of change of the dye mass is given by 


dM 
dt 


or rewriting in terms of u(x,t): 


— (x1, t) _ J (x2, t), 


ur(x,t)Ax = D (ux(Xo,t) — ux (x1, 6)). 
Dividing by Ax and taking the limit Ax — 0 results in the diffusion equation: 
ut = Duyy. 


We note that the diffusion equation is identical to the heat conduction equation, 
where u is temperature, and the constant D (commonly written as x) is the thermal 
conductivity. 


8.2 Derivation of the wave equation 


To derive the wave equation in one spacial dimension, we imagine an elastic string 
that undergoes small amplitude transverse vibrations. We define u(x,t) to be the 
vertical displacement of the string from the x-axis at position x and time t, and we 
wish to find the pde satisfied by u. We define p to be the constant mass density of 
the string, T the tension of the string, and @ the angle between the string and the 
horizonal line. We consider an infinitesimal string element located between x; and 
x2, with Ax = x2 — x1, as shown in Fig. 8.2. The governing equations are New- 
ton’s law of motion for the horizontal and vertical acceleration of our infinitesimal 
string element, and we assume that the string element only accelerates vertically. 
Therefore, the horizontal forces must balance and we have 


Ty cos 69 = T; cos 6. 


The vertical forces result in a vertical acceleration, and with uj; the vertical accel- 
eration of the string element and pV Ax? + Au2 = pAxy/1+ u2 its mass, where we 
have used uy = Au/Ax, exact as Ax + 0, we have 


pAxy/ 1+ uur = Th sin 6) — T; sin 6}. 
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Xx, X, 


Figure 8.2: Derivation of the wave equation. 


We now make the assumption of small vibrations, that is Au < Ax, or equivalently 
uy < 1. Note that [u] = / so that u; is dimensionless. With this approximation, to 
leading-order in uy we have 


cos 6) = cos 6; = 1, 


sin 02 = Ux (Xo, t), sin 6; = Ux (x1, t), 


V14+u2=1. 


Therefore, to leading order Tj = Tp = T, (ie., the tension in the string is approxi- 
mately constant), and 


and 


pAxupp = T (ux(X2,t) — ux (x1,t)). 
Dividing by Ax and taking the limit Ax — 0 results in the wave equation 
iC thee 


where c? = T/p. Since [T] = ml/t? and [p] = m/I1, we have [c*] = [7/#? so that c 
has units of velocity. 


8.3 Fourier series 
View tutorial on YouTube 


Our solution of the diffusion and wave equations will require use of a Fourier series. 
A periodic function f(x) with period 2L, can be represented as a Fourier series in 
the form 


ao. nTX _ nmx 
f(x) = 7 +L (ancos r by, sin L i (8.2) 
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Determination of the coefficients ag,41,a2,... and b1,b2,b3,... makes use of or- 
thogonality relations for sine and cosine. We first define the widely used Kronecker 


delta 64, as 
5 = 1 ifn=m,; 
mm") 0 otherwise. 


The orthogonality relations for n and m positive integers are then given with com- 
pact notation as the integration formulas 


[ : cos (™™*) cos (“*) dx = Loum, (8.3) 
iE sin (=) sin (“) dx = Léym, (8.4) 
ie ae (=) sin (WS) ax =f (8.5) 


To illustrate the integration technique used to obtain these results, we derive (8.4) 
assuming that n and m are positive integers with n #4 m. Changing variables to 
¢ = 7x/L, we obtain 


[0 En (2a 

=f" sin (mé) sin (n@)dé 

- =f. [cos ((m —n)&) —cos ((m+n)é)] de 
L 1 


TU 


: 1 ; 
= 5 [gee sin ((m — me) ~ SE sin ((m + m)2) 


—7 


ie sin? (“ax ~ 7 o ad are 


— 7 [ (1 — cos (2ng)) dg 


7 


~ On 2n 


=L. 


i E a= sin2ng| 7 


Integration formulas given by (8.3) and (8.5) can be similarly derived. 

To determine the coefficient a,, we multiply both sides of (8.2) by cos (n7tx/L) 
with n a nonnegative integer, and change the dummy summation variable from n 
to m. Integrating over x from —L to L and assuming that the integration can be 
done term by term in the infinite sum, we obtain 


L 7x AQ 
—dx = — —d 
f (x) cos L 5 4 
(oo) L L 
+ XL {n / cos _ cos ae dx + Diy cos sadlacs sin ae ax} 
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If n = 0, then the second and third integrals on the right-hand-side are zero and the 
first integral is 2L so that the right-hand-side becomes Lap. If n is a positive integer, 
then the first and third integrals on the right-hand-side are zero, and the second 
integral is Léym. For this case, we have 


[. f(x ) cos dx = = s Lamonm 


m=1 


= Lay, 


where all the terms in the summation except m = n are zero by virtue of the 
Kronecker delta. We therefore obtain for n = 0,1,2,... 


1 rt nNTXx 
ayn = L i f(x) cos a ie (8.6) 
To determine the coefficients b1, bz, b3,..., we multiply both sides of (8.2) by sin (n7tx/L), 


with n a positive integer, and again change the dummy summation variable from n 
to m. Integrating, we obtain 


eo f(x) sin a dx =Ff, ei aa dx 


L 
+ s {om [, sin ~ cos mnt dx + Din [sin = sin ma 


m=1 


Here, the first and second integrals on the right-hand-side are zero, and the third 
integral is Lénm so that 


a fi x) sin “ax = ‘3 LbinOnm 
m=1 


= Lby. 


Hence, for n = 1,2,3,..., 


=F = F(x )sin Tae. (8.7) 


Our results for the Fourier series of a function f(x) with period 2L are thus given 
by (8.2), (8.6) and (8.7). 


8.4 Fourier sine and cosine series 
View tutorial on YouTube 


The Fourier series simplifies if f(x) is an even function such that f(—x) = f(x), or 
an odd function such that f(—x) = —f(x). Use will be made of the following facts. 
The function cos (n7tx/L) is an even function and sin (n7tx/L) is an odd function. 
The product of two even functions is an even function. The product of two odd 
functions is an even function. The product of an even and an odd function is an 
odd function. And if g(x) is an even function, then 


[. g(x)dx = a g(x)dx 
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and if g(x) is an odd function, then 


[ seoae = 0. 


We examine in turn the Fourier series for an even or an odd function. First, if 
f(x) is even, then from (8.6) and (8.7) and our facts about even and odd functions, 


2 rh nTex 
ayn = ei, f(x) COs zh: 
bn = 0. 


(8.8) 


The Fourier series for an even function with period 2L is thus given by the Fourier 
cosine series 


m= NIX 
f(x) = at pE ay COS ——, f (x) even. (8.9) 
Second, if f(x) is odd, then 
an = 0, 
2. fF __ mmx (8.10) 
by ae f(x) sin ax; 


and the Fourier series for an odd function with period 2L is given by the Fourier 
sine series 


f(x) = 3 by sin ae f(x) odd. (8.11) 
n=1 


Examples of Fourier series computed numerically can be obtained using the 
Java applet found at http://www.falstad.com/fourier. Here, we demonstrate an 
analytical example. 


Example: Determine the Fourier cosine series of the even triangle function rep- 
resented by Fig. 8.3. 

View tutorial on YouTube 

The triangle function depicted in Fig. 8.3 is an even function of x with period 27r 
(i.e., L = 70). Its definition on 0 < x < 7 is given by 


f(x) =1- =. 


7U 


Because f(x) is even, it can be represented by the Fourier cosine series given by 
(8.8) and (8.9). The coefficient ag is 


a9 == (x)dx 
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Figure 8.3: The even triangle function. 


Note that a9 /2 is the average value of f(x) over one period, and it is obvious from 
Fig. 8.3 that this value is zero. The coefficients for n > 0 are 


An 


Since 


we have 


a = f" F(x) cos (nx)dx 


7 
= =a (1 = =) cos (nx)dx 
m Jo 1 


4 7 d 
al x cos (nx)dx 


2 {psme)], 


lI 
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cos (n7) —1, ifn odd; 
_ 1, if even; 
8/(n*77), if n odd; 
7) 0, if n even. 


. I sin (nx)ax} 
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The Fourier cosine series for the triangle function is therefore given by 


8 ' cos 3x cos 5x 
f(x) = a2 \cosx+ ay + ay te | 


Convergence of this series is rapid. As an interesting aside, evaluation of this series 
at x = 0, using f(0) = 1, yields an infinite series for m2 /8: 


With Fourier series now included in our applied mathematics toolbox, we are ready 
to solve the diffusion and wave equations in bounded domains. 


8.5 Solution of the diffusion equation 


8.5.1 Homogeneous boundary conditions 


We consider one dimensional diffusion in a pipe of length L, and solve the diffusion 
equation for the concentration u(x,t), 


up = Duyy, OS x<L, t>0. (8.12) 


Both initial and boundary conditions are required for a unique solution. That is, we 
assume the initial concentration distribution in the pipe is given by 


u(x,0) = f(x), O<x<L. (8.13) 


Furthermore, we assume that boundary conditions are given at the ends of the 
pipes. When the concentration value is specified at the boundaries, the boundary 
conditions are called Dirichlet boundary conditions. As the simplest example, we as- 
sume here homogeneous Dirichlet boundary conditions, that is zero concentration 
of dye at the ends of the pipe, which could occur if the ends of the pipe open up 
into large reservoirs of clear solution, 


u(0,t)=0, u(L,t)=0, t>0. (8.14) 


We will later also discuss inhomogeneous Dirichlet boundary conditions and homo- 
geneous Neumann boundary conditions, for which the derivative of the concentration 
is specified to be zero at the boundaries. Note that if f(x) is identically zero, then 
the trivial solution u(x,t) = 0 satisfies the differential equation and the initial and 
boundary conditions and is therefore the unique solution of the problem. In what 
follows, we will assume that f(x) is not identically zero so that we need to find a 
solution different than the trivial solution. 

The solution method we use is called separation of variables. We assume that 
u(x,t) can be written as a product of two other functions, one dependent only on 
position x and the other dependent only on time t. That is, we make the ansatz 


u(x,t) = X(x)T(t). (8.15) 


Whether this ansatz will succeed depends on whether the solution indeed has this 
form. Substituting (8.15) into (8.12), we obtain 


XT’ = DX"T, 
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which we rewrite by separating the x and ¢ dependence to opposite sides of the 
equation: 
x" 1 T’ 
xe oF 
The left hand side of this equation is independent of t and the right hand side is 
independent of x. Both sides of this equation are therefore independent of both x 
and ¢ and equal to a constant. Introducing —A as the separation constant, we have 
xl 1 T! 
— SS —A 
Xx DT , 
and we obtain the two ordinary differential equations 
XxX” +AX=0, T'+ADT=0. (8.16) 


Because of the boundary conditions, we must first consider the equation for X(x). 
To solve, we need to determine the boundary conditions at x = 0 and x = L. Now, 
from (8.14) and (8.15), 


u(0,t) = X(0)T(t)=0, t>0. 


Since T(t) is not identically zero for all t (which would result in the trivial solution 
for u), we must have X(0) = 0. Similarly, the boundary condition at x = L requires 
X(L) = 0. We therefore consider the two-point boundary value problem 


X”+AX=0, X(0) = X(L) =0. (8.17) 


The equation given by (8.17) is called an ode eigenvalue problem. Clearly, the 
trivial solution X(x) = 0 is a solution. Nontrivial solutions exist only for discrete 
values of A. These discrete values of A and the corresponding functions X(x) are 
called the eigenvalues and eigenfunctions of the differential equation. 

Since the form of the general solution of the ode depends on the sign of A, we 
consider in turn the cases A > 0, A < 0 and A = 0. For A > 0, we write A = we and 
determine the general solution of 


X" +X =0 
to be 
X(x) = Acos px + Bsin px. 


Applying the boundary condition at x = 0, we find A = 0. The boundary condition 
at x = L then yields 
Bsin wl = 0. 


The solution B = 0 results in the trivial solution for u and can be ruled out. There- 
fore, we must have 
sinul =0, 


which is an equation for jy. The solutions are 
w=nr/L, 
where n is an integer. We have thus determined the eigenvalues A = we > 0 to be 


An =(nn/LY, n=1,2,3,..., (8.18) 
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with corresponding eigenfunctions 
Xy = sin (n7x/L). (8.19) 
For A < 0, we write A = —v and determine the general solution of 
xX" — WX =0 
to be 
X(x) = Acoshyx + Bsinh px, 


where we have previously introduced the hyperbolic sine and cosine functions in 
§3.4.1. Applying the boundary condition at x = 0, we find A = 0. The boundary 
condition at x = L then yields 


Bsinh pL = 0, 


which for 7» # 0 has only the solution B = 0. Therefore, there is no nontrivial 
solution for u with A < 0. Finally, for A = 0, we have 


xl = 0, 


with general solution 
X(x) = A+ Bx. 


The boundary condition at x = 0 and x = L yields A = B = 0 so again there is no 
nontrivial solution for uv with A = 0. 

We now turn to the equation for T(t). The equation corresponding to the eigen- 
value A;,, using (8.18), is given by 


7! (n?n?D/L?) T=0, 
which has solution proportional to 
Ty = eo DEL (8.20) 


Therefore, multiplying the solutions given by (8.19) and (8.20), we conclude that the 
functions . 

n(x, t) = sin (n7tx/L)e~* 7 Pt/E (8.21) 
satisfy the pde given by (8.12) and the boundary conditions given by (8.14) for every 
positive integer 7. 

The principle of linear superposition for homogeneous linear differential equa- 
tions then states that the general solution to (8.12) and (8.14) is given by 


u(x,t) = y byun(x, t) 
sa (8.22) 
-_ : —n?2n?Dt/L? 
= VY" by sin (n7x/L)e : 


n=1 


The final solution step is to satisfy the initial conditions given by (8.13). At t = 0, 
we have 7 
f(x) = YO by sin (n7tx/L). (8.23) 
n 


=1 
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We immediately recognize (8.23) as a Fourier sine series (8.11) for an odd function 
f(x) with period 2L. Equation (8.23) is a periodic extension of our original f(x) 
defined on 0 < x < L, and is an odd function because of the boundary condition 
f(0) = 0. From our solution for the coefficients of a Fourier sine series (8.10), we 
determine 


2 rh _ n7x 
by = = [ f(x) sin a. (8.24) 


Thus the solution to the diffusion equation with homogeneous Dirichlet boundary 
conditions defined by (8.12), (8.13) and (8.14) is given by (8.22) with the b;, coeffi- 
cients computed from (8.24). 


Example: Determine the concentration of a dye in a pipe of length L, where the 
dye has total initial mass Mo and is initially concentrated at the center of the 
pipe, and the ends of the pipe are held at zero concentration. 


The governing equation for concentration is the diffusion equation. We model the 
initial concentration of the dye by a delta-function centered at x = L/2, that is, 
u(x,0) = f(x) = Mod(x — L/2). Therefore, from (8.24), 


Doig L, . nmx 
by = = Mod(x — 5) sin ——dx 


= a sin (n7/2) 


2Mo/L ifn =1,5,9,...; 
= ¢ —2M)/L ifn =3,7,11,...; 
0 ifn =2,4,6,.... 


With b,, determined, the solution for u(x,t) given by (8.22) can be written as 


2 ee 2 1 
u(x,t) = a y\(-1)"sin (fe) eo 2n+1)272Dt/12_ 
n=0 


When t >> L?/D, the leading-order term in the series is a good approximation and 
is given by 


0 


2 
u(x,t) © a sin (7x/L)e~™ Pt/E* 


The mass of the dye in the pipe is decreasing in time, diffusing into the reservoirs 
located at both ends. The total mass in the pipe at time f can be found from 


M(t) = a u(x, t)dx, 
and since 
a sin (77x/L)dx = ad 
0 aa 
we have for large times 
M(t) = M0 om DH/1?, 


7U 
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8.5.2 Inhomogeneous boundary conditions 


Consider a diffusion problem where one end of the pipe has dye of concentration 
held constant at C; and the other held constant at C2, which could occur if the ends 
of the pipe had large reservoirs of fluid with different concentrations of dye. With 
u(x,t) the concentration of dye, the boundary conditions are given by 


u(0,t)=Cy, u(L,t)=Co, t>0. 
The concentration u(x,t) satisfies the diffusion equation with diffusivity D: 
ut = Duyy. 


If we try to solve this problem directly using separation of variables, we will run 
into trouble. Applying the inhomogeneous boundary condition at x = 0 directly to 
the ansatz u(x,t) = X(x)T(t) results in 


u(0,t) = X(0)T(t) = Cy; 


so that 

X(0) = C,/T(t). 
However, our separation of variables ansatz assumes X(x) to be independent of t! 
We therefore say that inhomogeneous boundary conditions are not separable. 

The proper way to solve a problem with inhomogeneous boundary conditions 
is to transform it into another problem with homogeneous boundary conditions. 
As t — oo, we assume that a stationary concentration distribution v(x) will attain, 
independent of t. Since v(x) must satisfy the diffusion equation, we have 

o'(x)=0, O<x<L, 
with general solution 
v(x) = A+ Bx. 
Since v(x) must satisfy the same boundary conditions of u(x,t), we have v(0) = C; 
and v(L) = C2, and we determine A = C; and B = (C) —C,)/L. 

We now express u(x,t) as the sum of the known asymptotic stationary con- 
centration distribution v(x) and an unknown transient concentration distribution 
w(x,t): 

u(x,t) = v(x) + w(x,t). 
Substituting into the diffusion equation, we obtain 


7) a 
2 (o(x) + w(x,#)) = D2 (ols) + 0(x,8)) 
Wt = DW yx, 


since vy; = 0 and vx; = 0. The boundary conditions satisfied by w are 

w(0,t) = u(0,t) — v(0) =0, 

w(L,t) = u(L,t) — o(L) =0, 
so that w is observed to satisfy homogeneous boundary conditions. If the initial 
conditions are given by u(x,0) = f(x), then the initial conditions for w are 

w(x,0) = u(x,0) — v(x) 
= f(x) — o(2). 

The resulting equations may then be solved for w(x,t) using the technique for ho- 
mogeneous boundary conditions, and u(x,t) subsequently determined. 
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8.5.3 Pipe with closed ends 


There is no diffusion of dye through the ends of a sealed pipe. Accordingly, the 
mass flux of dye through the pipe ends, given by (8.1), is zero so that the boundary 
conditions on the dye concentration u(x,t) becomes 


ux(0,t)=0, ux(L,t)=0, t>0, (8.25) 


which are known as homogeneous Neumann boundary conditions. Again, we ap- 
ply the method of separation of variables and as before, we obtain the two ordinary 
differential equations given by (8.16). Considering first the equation for X(x), the 
appropriate boundary conditions are now on the first derivative of X(x), and we 
must solve 

X"+AK=0, X'(0)=X'(L)=0. (8.26) 


Again, we consider in turn the cases A > 0, A < 0 and A = 0. For A > 0, we write 
A= we and determine the general solution of (8.26) to be 


X(x) = Acos px + Bsin px, 
so that taking the derivative 
X' (x) = —pAsin px + vB cos px. 


Applying the boundary condition X’(0) = 0, we find B = 0. The boundary condi- 
tion at x = L then yields 
—pAsinuL = 0. 


The solution A = 0 results in the trivial solution for u and can be ruled out. There- 
fore, we must have 
sinuL =0, 


with solutions 
pw=nr/L, 


where nis an integer. We have thus determined the eigenvalues A = ie > 0 to be 
An =(nn/L)*, n=1,2,3,..., (8.27) 
with corresponding eigenfunctions 
Xn = cos (n7x/L). (8.28) 
For A < 0, we write A = —y and determine the general solution of (8.26) to be 
X(x) = Acoshyx + Bsinh px, 


so that taking the derivative 


X’(x) = pA sinh px + B cosh px. 


Applying the boundary condition X’(0) = 0 yields B = 0. The boundary condition 
X'(L) = 0 then yields 
vA sinh pL = 0, 


which for pt: # 0 has only the solution A = 0. Therefore, there is no nontrivial 
solution for u with A < 0. Finally, for A = 0, the general solution of (8.26) is 


X(x) = A+ Bx, 
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so that taking the derivative 
X'(x) = B. 


The boundary condition X’(0) = 0 yields B = 0; X’(L) = 0 is then trivially satisfied. 
Therefore, we have an additional eigenvalue and eigenfunction given by 


Ag =0, Xo(x) =1, 


which can be seen as extending the formula obtained for eigenvalues and eigenvec- 
tors for positive A given by (8.27) and (8.28) to n = 0. 

We now turn to the equation for T(t). The equation corresponding to eigenvalue 
An, using (8.27), is given by 


ae (n?n?D/L?) P20, 
which has solution proportional to 
Ty = eT DEL? (8.29) 


valid for n = 0,1,2,.... Therefore, multiplying the solutions given by (8.28) and 
(8.29), we conclude that the functions 


Un(x,t) = cos (nzx/L)e™ ™ Dt/L* (8.30) 


satisfy the pde given by (8.12) and the boundary conditions given by (8.25) for every 
nonnegative integer n. 
The principle of linear superposition then yields the general solution as 


u(x,t) =} tun (%,t) 
“ ~ vs (8.31) 
= 4 + y an cos (n7x/L)e~" ™ cae 
n=1 


where we have redefined the constants so that cp = ag/2 and cy = ay,n = 1,2,3,.... 
The final solution step is to satisfy the initial conditions given by (8.13). At t = 0, 
we have 


lee} 
f(x) = * + Y° an cos (n7tx/L), (8.32) 
n=1 
which we recognize as a Fourier cosine series (8.9) for an even function f(x) with 
period 2L. We have obtained a cosine series for the periodic extension of f(x) 
because of the boundary condition f'(0) = 0, which is satisfied by an even function 
with continuous first derivative. From our solution (8.8) for the coefficients of a 
Fourier cosine series, we determine 


2 rh nTex 
an = = [ f(x) cos dx, (8.33) 


Thus the solution to the diffusion equation with homogenous Neumann boundary 
conditions defined by (8.12), (8.13) and (8.25) is given by (8.31) with the coefficients 
computed from (8.33). 
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Example: Determine the concentration of a dye in a pipe of length L, where the 
dye has total initial mass Mo and is initially concentrated at the center of the 
pipe, and the ends of the pipe are sealed 


Again we model the initial concentration of the dye by a delta-function centered at 
x = L/2. From (8.33), 


— ah Mod (x — 5) cos ay 


2 
= a cos (n7/2) 


2Mo/L ifn =0,4,8,...; 
= 4 —2Mo/L ifn=2,6,10,...; 
0 ifn = 1,3,5,.... 


The first two terms in the series for u(x,t) are given by 
u(x,t) = a fa —2cos (27x /L)e 4 Dt/L? +... | : 


Notice that as tf + oo, u(x,t) + Mg/L: the dye mass is conserved in the pipe (since 
the pipe ends are sealed) and eventually diffused uniformly throughout the pipe of 
length L. 


8.6 Solution of the wave equation 


8.6.1 Plucked string 


We assume an elastic string with fixed ends is plucked like a guitar string. The gov- 
erning equation for u(x,t), the position of the string from its equilibrium position, 
is the wave equation 

Uy = Czy, (8.34) 


with c2 = T/ p and with boundary conditions at the string ends located at x = 0 
and L given by 

u(0,t) =0, u(L,t) =0. (8.35) 
Since the wave equation is second-order in time, initial conditions are required for 
both the displacement of the string due to the plucking and the initial velocity of 
the displacement. We assume 


u(x,0) = f(x), u(x,0)=0, O<x<L. (8.36) 

Again we use the method of separation of variables and try the ansatz 
u(x,t) = X(x)T(t). (8.37) 
Substitution of our ansatz (8.37) into the wave equation (8.34) and separating vari- 


ables results in 
x" 1 qT" 


xX Tapa 
yielding the two ordinary differential equations 
X"4+AX=0, T’+Ac?T =0. (8.38) 
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We solve first the equation for X(x). The appropriate boundary conditions for X 
are given by 
X(0)=0, X(L)=0, (8.39) 


and we have solved this equation for X(x) previously in §8.5.1 (see (8.17)). A non- 
trivial solution exists only when A > 0, and our previously determined solution 
was 

An =(nn/L)Y, n=1,2,3,..., (8.40) 


with corresponding eigenfunctions 
Xn = sin (n7x/L). (8.41) 


With A, specified, the T equation then becomes 


n212¢2 
ae =r 12 Th = 0, 
with general solution given by 
t t 
Tr(t) = Acos = + Bsin = : (8.42) 


The second of the initial conditions given by (8.36) implies 
u;(x,0) = X(x)T’(0) =0, 


which can be satisfied only if T’(0) = 0. Applying this boundary condition to 
(8.42), we find B = 0. Combining our solution for X;,(x), (8.41), and T,,(t), we have 
determined that 

NTUX n7tct 


Un(x,t) = sin L cos Lt n=1,2,3,... 


satisfies the wave equation, the boundary conditions at the string ends, and the 
assumption of zero initial string velocity. Linear superposition of these solutions 
results in the general solution for u(x,t) of the form 


kes t 
u(x,t) = )7 by sin — cos — : (8.43) 


n=1 


The remaining condition to satisfy is the initial displacement of the string, the first 
equation of (8.36). We have 


CO 


f(x) = YC by sin (n7x/L), 


n=1 


which is observed to be a Fourier sine series (8.11) for an odd function with period 
2L. Therefore, the coefficients b, are given by (8.10), 


L 
= al f(x)sin “dx, 1 =1,2,3,.... (8.44) 
Our solution to the wave equation with plucked string is thus given by (8.43) and 
(8.44). Notice that the solution is time periodic with period 2L/c. The correspond- 
ing fundamental frequency is the reciprocal of the period and is given by f = c/2L. 


From our derivation of the wave equation in §8.2, the velocity c is related to the 
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density of the string p and tension of the string T by c= T/ p. Therefore, the 
fundamental frequency (pitch) of our “guitar string” increases (is raised) with in- 
creasing tension, decreasing string density, and decreasing string length. Indeed, 
these are exactly the parameters used to construct, tune and play a guitar. 

The wave nature of our solution and the physical significance of the velocity c 
can be made more transparent if we make use of the trigonometric identity 


1 
sin x cosy = 5 (sin (x + y) +sin(x—y)). 
With this identity, our solution (8.43) can be rewritten as 


(x,t) = a by (sin wn a3, + sin wee *) : (8.45) 


The first and second sine functions can be interpreted as a traveling wave moving 
to the left or the right with velocity c. This can be seen by incrementing time, 
t — t+, and observing that the value of the first sine function is unchanged 
provided the position is shifted by x — x — cd, and the second sine function is 
unchanged provided x + x + cd. Two waves travelling in opposite directions with 
equal amplitude results in a standing wave. 


8.6.2 Hammered string 


In contrast to a guitar string that is plucked, a piano string is hammered. The 
appropriate initial conditions for a piano string would be 


u(x,0)=0, u(x,0) =g9(x), O<x <L. (8.46) 


Our solution proceeds as previously, except that now the homogeneous initial con- 
dition on T(t) is T(0) = 0, so that A = 0 in (8.42). The wave equation solution is 


therefore 


foe) 
t 
u(x,t) = J bn sin ee sin = (8.47) 
n=1 


Imposition of initial conditions then yields 


NTTX 


g(x) = — = nb, sin a 
The coefficient of the Fourier sine series for g(x) is seen to be n7tcb,/L, and we have 


ia 
a =? fat ) sin we dx, 


or 


8.6.3 General initial conditions 


If the initial conditions on u(x,t) are generalized to 
u(x,0) = f(x), m(x,0)=g(x), OS¥<L, (8.48) 
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Figure 8.4: Dirichlet problem for the Laplace equation in a rectangle. 


then the solution to the wave equation can be determined using the principle of 
linear superposition. Suppose v(x,t) is the solution to the wave equation with 
initial condition (8.36) and w(x,t) is the solution to the wave equation with initial 
conditions (8.46). Then we have 


u(x,t) = v(x,t) + w(x,t), 


since u(x,t) satisfies the wave equation, the boundary conditions, and the initial 
conditions given by (8.48). 


8.7 The Laplace equation 
The diffusion equation in two spatial dimensions is 


up = D(uxx + Uyy). 


The steady-state solution, approached asymptotically in time, has u; = 0 so that the 
steady-state solution u = u(x,y) satisfies the two-dimensional Laplace equation 


Uxx + Uyy = 0. (8.49) 


We will consider the mathematical problem of solving the two dimensional Laplace 
equation inside a rectangular or a circular boundary. The value of u(x,y) will be 
specified on the boundaries, defining this problem to be of Dirichlet type. 


8.7.1 Dirichlet problem for a rectangle 


We consider the Laplace equation (8.49) for the interior of a rectangle 0 < x <a, 
0<y <b, (see Fig. 8.4), with boundary conditions 


u(x,0) =0, u(x,b) =0, 0<x<@; 
u(0,y) =0, u(a,y) = fly), O<y<b. 
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More general boundary conditions can be solved by linear superposition of solu- 
tions. 
We take our usual ansatz 


u(x,y) = X(x)Y(y), 


and find after substitution into (8.49), 


with A the separation constant. We thus obtain the two ordinary differential equa- 
tions 


X"-AX=0, Y"+AY=0. 


The homogeneous boundary conditions are X(0) = 0, Y(0) = 0 and Y(b) = 0. 
We have already solved the equation for Y(y) in §8.5.1, and the solution yields the 
eigenvalues 


= (")’, a ec 


with corresponding eigenfunctions 


Yn(y) = sin ae 


The remaining X equation and homogeneous boundary condition is therefore 


n? 72 
Xf ga X(0) =0, 
and the solution is the hyperbolic sine function 


X,(x) = sinh =, 


times a constant. Writing u;, = XnYn, multiplying by a constant and summing over 
n, yields the general solution 


u(x,y) = Lon sinh “a sin a 


The remaining inhomogeneous boundary condition u(a,y) = f(y) results in 


=F cy sinh gin 7Y 
Fy) )_ Greil 7 sine, 


which we recognize as a Fourier sine series for an odd function with period 2b, and 
coefficient c, sinh (n7ta/b). The solution for the coefficient is given by 


2 b __ n7y 
oo Fame fy usin “dy. 
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8.7.2 Dirichlet problem for a circle 


The Laplace equation is commonly written symbolically as 
V2u =0, (8.50) 


where V7? is called the Laplacian, sometimes denoted as A. The Laplacian can be 
written in various coordinate systems, and the choice of coordinate systems usually 
depends on the geometry of the boundaries. Indeed, the Laplace equation is known 
to be separable in 13 different coordinate systems! We have solved the Laplace 
equation in two dimensions, with boundary conditions specified on a rectangle, 
using 

a a 

ax2 t oy? 

Here we consider boundary conditions specified on a circle, and write the Lapla- 
cian in polar coordinates by changing variables from cartesian coordinates. Polar 
coordinates is defined by the transformation (r,6) — (x,y): 


oS 


x=rcosé, y=rsiné; (8.51) 
and the chain rule gives for the partial derivatives 


du ouodx  dudy Ou _ dudx , dudy 
or oxor  oydar’ 00 dxd0 | dyad 


(8.52) 


After taking the partial derivatives of x and y using (8.51), we can write the trans- 
formation (8.52) in matrix form as 


du/or \ _ cos @ sin 0 ou/ox 
( ou/0e ) — ( —rsin@ rcosé ) ( ou/oy ji (8.53) 


Inversion of (8.53) can be determined from the following result, commonly proved 
in a linear algebra class. If 


a b 
a=(% a detA #0, 


1 d —b 
=] _ 
= ere ( =€ ) 


Therefore, since the determinant of the 2 x 2 matrix in (8.53) is r, we have 


then 


du/dx \ _ ( cos@ —sin@/r ou/or (8.54) 
du/dy /  \ sin@ cos 0/r ou/oe )- . 
Rewriting (8.54) in operator form, we have 
) d sind d ) . 79 , cosé d 
5x = 00895, , oO (Oy = sind + —— ap. (8.55) 


To find the Laplacian in polar coordinates with minimum algebra, we combine 
(8.55) using complex variables as 


z rig =e8(S 425), (8.56) 
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so that the Laplacian may be found by multiplying both sides of (8.56) by its com- 
plex conjugate, taking care with the computation of the derivatives on the right- 


hand-side: 
a2 a2 ao (9  i0\ »#f9d iad 
=¢ boas be” | aa on 
ax? oy? or roe or roe 


oe 1 eo Te 
7 36? 


ore | ror. 
We have therefore determined that the Laplacian in polar coordinates is given by 


2 10 10 
(8.57) 


Ze 
¥ ~ or2 ' ror | 12062’ 


which is sometimes written as 


2 
ve 10 y re) | 10 
ror \ or 12 062 
We now consider the solution of the Laplace equation in a circle with radius 
r <a subject to the boundary condition 


u(a,0) = f(0), O<@< 27. (8.58) 


An additional boundary condition due to the use of polar coordinates is that u(r, @) 
is periodic in @ with period 27. Furthermore, we will also assume that u(r,@) is 
finite within the circle. 

The method of separation of variables takes as our ansatz 


u(r,0) = R(r)@(8), 


and substitution into the Laplace equation (8.50) using (8.57) yields 


1 1 
R"© + —R’O4+ —=RO” =0, 
r 72 


a " / /} 
Fas | Ps = S =A, 
R R © 
where A is the separation constant. We thus obtain the two ordinary differential 


equations 


r7>R” +7R'-AR=0, ©” +AO=0. 


The © equation is solved assuming periodic boundary conditions with period 27:. 
If A < 0, then no periodic solution exists. If A = 0, then © can be constant. If 
A= e > 0, then 

©(6) = Acos pO + Bsin p60, 


and the requirement that © is periodic with period 27 forces y to be an integer. 
Therefore, 
An=n’*, n=0,1,2,..., 


with corresponding eigenfunctions 


©,(0) = Ancosné + B, sinné. 
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The R equation for each eigenvalue A, then becomes 
77R" +7R! —n?R=0, (8.59) 


which is an Euler equation. With the ansatz R = r°, (8.59) reduces to the algebraic 
equation s(s — 1) +s — n* = 0, or s* = n?. Therefore, s = +n, and there are two 
real solutions when n > 0 and degenerate solutions when n = 0. When n > 0, the 
solution for R(r) is 


Ry(r) = Ar" + Bro”. 


The requirement that u(r,@) is finite in the circle forces B = 0 since r~” becomes 
unbounded as r + 0. When n = 0, the solution for R(r) is 


Ri(r) =A+Blnr, 


and again finite u in the circle forces B = 0. Therefore, the solution for n = 0,1,2,... 
is R, =r”. Thus the general solution for u(r,@) may be written as 


A foe) 
u(r,0) = - + )7 r"(Ancosné + By sinn6), (8.60) 
n=1 


where we have separated out the n = 0 solution to write our solution in a form 
similar to the standard Fourier series given by (8.2). The remaining boundary con- 
dition (8.58) specifies the values of u on the circle of radius a, and imposition of this 
boundary condition results in 


f(@) = at + > a" (Ay cosné + By sinné). (8.61) 
n=1 


Equation (8.61) is a Fourier series for the periodic function f(@) with period 27, ie., 
L = 7 in (8.2). The Fourier coefficients a" A, and a"B, are therefore given by (8.6) 
and (8.7) to be 


a" An 


© [" Fo) cosmpag, n= 0,1, 2g 


a" Bn 


7 i wo = =e (8.62) 


where we have used ¢ for the dummy variable of integration. 
A remarkable fact is that the infinite series solution for u(r,@) can be summed 
explicitly. Substituting (8.62) into (8.60), we obtain 


aN 


1 27 co r\n . . 
u(r,0) = at dpf(¢) f +2)° “) cost cosn + inne sinnp) 


=x | aor) f +2, 


i. 


“\" cos n(@ — ’) : 


We can sum the infinite series by writing 2. cos n(@ — p) = ei(°—-#) + e~m(0—-9), and 
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using the sum of the geometric series )(_, 2" = z/(1—) to obtain 
OO rpyn 0 f rei(0-9)\" 0 pe i(0-9) \" 
142)" (-) cosn(@—) =1 a oe 
n=1 a a a 


1 rei(9—-) 
= i= reil0—$) rT C.C. 


a — 7? 


~ a2 — 2arcos (0 — ) + 12" 


Therefore, 


a i a f(p) 
ane) Qn i a* — 2ar cos (@ — p) +r? ae 


an integral result for u(r,@) known as Poisson’s formula. As a trivial example, 
consider the solution for u(r,0) if f(@) = F, a constant. Clearly, u(r,@) = F satisfies 
both the Laplace equation and the boundary conditions so must be the solution. 
You can verify that u(r,@) = F is indeed the solution by showing that 


a dp ee: 
0 


a2 —2arcos(@—@p) +r? 9 at — 1? 


8.8 The Schrédinger equation 


The Schrédinger equation is the equation of motion for nonrelativistic quantum 
mechanics. This equation is a linear partial differential equation and in simple 
situations can be solved using the technique of separation of variables. Luckily, one 
of the cases that can be solved analytically is the hydrogen atom. It can be argued 
that the greatest success of classical mechanics was the solution of the Earth-Sun 
system. Similarly, it can be argued that the greatest success of quantum mechanics 
was the solution of the electron-proton system. The mathematical solutions of these 
classical and quantum two-body problems rank among the highest achievements of 
mankind. 


8.8.1 Heuristic derivation of the Schrédinger equation 


Nature consists of waves and particles. Early in the twentieth century, light was 
discovered to act as both a wave and a particle, and quantum mechanics assumes 
that matter can also act in both ways. 

In general, waves are described by their wavelength A and frequency v, or equiv- 
alently their wavenumber k = 27t/A and their angular frequency w = 27v. The 
Planck-Einstein relations for light postulated that the energy E and momentum p 
of a light particle, called a photon, is proportional to the angular frequency w and 
wavenumber k of the light wave. The proportionality constant is called “h bar”, 
denoted as fi, and is related to the original Planck’s constant h by h = h/27. The 
Planck-Einstein relations are given by 


E=hw, p=hk. 
De Broglie in 1924 postulated that matter also follows these relations. 
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A heuristic derivation of the Schrodinger equation for a particle of mass m and 
momentum p constrained to move in one dimension begins with the classical equa- 


tion 
p 

+ 4V(x,t) =E, 8.63 
P V(x) (8.63) 
where p”/2m is the kinetic energy of the mass, V(x,t) is the potential energy, and 
E is the total energy. 

In search of a wave equation, we consider how to write a free wave in one 
dimension. Using a real function, we could write 


Y = Acos (kx -—wt+¢), (8.64) 


where A is the amplitude and @ is the phase. Or using a complex function, we 
could write 

W = Cellkx-wt) (8.65) 
where C is a complex number containing both amplitude and phase. 


We now rewrite the classical energy equation (8.63) using the Planck-Einstein 
relations. After multiplying by a wavefunction, we have 


7 
ak E(x, th) + V(x, t)¥ (x,t) = hw¥ (x,t). 


We would like to replace k and w, which refer to the wave characteristics of the 
particle, by differential operators acting on the wavefunction ¥ (x,t). If we consider 
V(x,t) = 0 and the free particle wavefunctions given by (8.64) and (8.65), it is easy 
to see that to replace both k? and w by derivatives we need to use the complex form 
of the wavefunction and explicitly introduce the imaginary unit i, that is 


a me) 
Ox?’ ot 


Doing this, we obtain the intrinsically complex equation 


ro 


hr ey ay 
—-— > 4+ V(x,t)¥ = ih— 

2m ox? Ve) eee 
which is the one-dimensional Schrédinger equation for a particle of mass m in a 
potential V = V(x,t). This equation is easily generalized to three dimensions and 
takes the form 


ie oy OF (x,t) 
—- aa Y (x,t) + V(x, t)¥ (x,t) = ih a” (8.66) 
where in Cartesian coordinates the Laplacian V? is written as 
v2 _ ay ay Ay 
Ox2 © dy2 | dz?" 
The Born interpretation of the wavefunction states that |¥(x, t)|? is the probabil- 


ity density function of the particle’s location. That is, the spatial integral of |‘¥ (x, f) |? 
over a volume V gives the probability of finding the particle in V at time ft. Since 
the particle must be somewhere, the wavefunction for a bound particle is usually 


normalized so that 
i. / / |¥ (x,y,z; t)|?dx dy dz = 1. 
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The simple requirements that the wavefunction be normalizable as well as single 
valued admits an analytical solution of the Schrodinger equation for the hydrogen 
atom. 


8.8.2 The time-independent Schrédinger equation 


The space and time variables of the time-dependent Schrédinger equation (8.66) 
can be separated provided the potential function V(x,t) = V(x) is independent of 
time. We try ¥ (x,t) = p(x) f(t) and obtain 


hn , 
—5— fVp + V(x)pf = ihpf. 
2m 
Dividing by wf, the equation separates as 


hn V7 = a 
mp =e 


The left-hand side is independent of t, the right-hand side is independent of x, so 
both the left- and right-hand sides must be independent of x and ¢ and equal to a 
constant. We call this separation constant E, thereby reintroducing the total energy 
into the equation. We now have the two differential equations 


, iE h 
7 af D) 


2 
Vp + V(x)p = Ey. 
m 
The second equation is called the time-independent Schrédinger equation. The first 
equation can be easily integrated to obtain 


f(t) = etEt/h 


which can be multiplied by a arbitrary constant. 


8.8.3 Particle in a one-dimensional box 


We assume that a particle of mass m is able to move freely in only one dimension 
and is confined to the region defined by 0 < x < L. This is perhaps the sim- 
plest quantum mechanical problem with quantized energy levels. We take as the 
potential energy function 


0, O L, 
V(x) = Salk 
co, otherwise, 
where we may assume that the particle is forbidden from the region with infinite 
potential energy. We will simply take as the boundary conditions on the wavefunc- 
tion 


p(0) = p(L) =0. (8.67) 


For this potential, the time-independent Schrédinger equation for 0 < x < L 
reduces to 
2 72 
Pep ey 
2m dx? 


CHAPTER 8. PARTIAL DIFFERENTIAL EQUATIONS 127 


8.8. THE SCHRODINGER EQUATION 


which we write in a familiar form as 


2 
ms ! (= ) p=0. (8.68) 


Equation (8.68) together with the boundary conditions (8.67) form an ode eigen- 
value problem, which is in fact identical to the problem we solved for the diffusion 
equation in subsection 8.5.1. 

The general solution to this second-order differential equation is given by 


V2mE V2mE 
(x) = Acos - * + Bsin - as 


The first boundary condition (0) = 0 yields A = 0. The second boundary condi- 


tion #(L) = 0 yields 
J2mEL _ 
7S 


The energy levels of the particle are therefore quantized, and the allowed values are 
given by 


na, n=1,2,3,.... 


nnn 
2mL2 * 


The corresponding wavefunction is given by p, = B sin (n7tx/L). We can normalize 
each wavefunction so that 


En — 


Ts 
| [Pn (x) |Pdx =1, 
obtaining B = /2/L. We have therefore obtained 


2 n7ex 


a(x) = psin"*, O<x<L,; 
i = 
0, otherwise. 


8.8.4 The simple harmonic oscillator 


Hooke’s law for a mass on a spring is given by 
F=-—kKx, 


where K is the spring constant. The potential energy V(x) in classical mechanics 
satisfies 
F = —0V/ox, 


so that the potential energy of the spring is given by 


1 
V(x) = =Kx?. 
(x) = 5Kx 
Recall that the differential equation for a classical mass on a spring is given from 
Newton’s law by 
mx = —Kx, 


which can be rewritten as 
H+ wx =0, 
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where w? = K/m. Following standard notation, we will therefore write the poten- 
tial energy as 
1 
V(x) = smu x. 
The time-independent Schrédinger equation for the one-dimensional simple har- 
monic oscillator then becomes 
wep 1 o> 


The relevant boundary conditions are p — 0 as x —> +co so that the wavefunction 
is normalizable. 

The Schrédinger equation given by (8.69) can be made neater by nondimension- 
alization. We rewrite (8.69) as 


ap 2mE  m*w*x? 

Ant ( v; 72 ) yp=0, (8.70) 
and observe that the dimension [h?/m?w*] = I*, where | is a unit of length. We 
therefore let 

h 
PN a p(x) = uly), 
and (8.70) becomes 
a ee (8.71) 
dy? at : 
with the dimensionless energy given by 
E =2E/hw, (8.72) 
and boundary conditions 
lim u(y) =0. (8.73) 
y 10 


The dimensionless Schrédinger equation given by (8.71) together with the bound- 
ary conditions (8.73) forms another ode eigenvalue problem. A nontrivial solution 
for u = u(y) exists only for discrete values of €, resulting in the quantization of the 
energy levels. Since the second-order ode given by (8.71) has a nonconstant coeffi- 
cient, we can use the techniques of Chapter 5 to find a convergent series solution 
for u = u(y) that depends on €. However, we will then be faced with the diffi- 
cult problem of determining the values of € for which u(y) satisfies the boundary 
conditions (8.73). 

A path to an analytical solution can be discovered if we first consider the behav- 
ior of u for large values of y. With y? >> €, (8.71) reduces to 

au 


re oe (8.74) 


To determine the behavior of u for large y, we try the ansatz 


2 


u(y) =e". 


We have ‘ 
u(y) = Zaye", 


CHAPTER 8. PARTIAL DIFFERENTIAL EQUATIONS 129 


8.8. THE SCHRODINGER EQUATION 


u"(y) =e (2a + 40°y?) we day, 
and substitution into (8.74) results in 


(4a? — 1)y* = 0, 


yielding a = +1/2. Therefore at large y, u(y) either grows like ev /2 of decays like 


ey /2, Here, we have neglected a possible polynomial factor in front of the ex- 
ponential functions. Clearly, the boundary conditions forbid the growing behavior 
and only allow the decaying behavior. 

We proceed further by letting 


u(y) = Hye? ?, 


and determining the differential equation for H(y). After some simple calculation, 
we have 


wy) = (Hi —yHeY?, uly) = (HY —2yH! + — A) EF”. 


Substitution of the second derivative and the function into (8.71) results in the dif- 
ferential equation 
H" —2yH'+(E€-1)H=0. (8.75) 


We now solve (8.75) by a power-series ansatz. We try 
H(y) = Yo ary’. 
k=0 


Substitution into (8.75) and shifting indices as detailed in Chapter 5 results in 


foe) 


Y" ((kK +2) (K+ 1)aep2 + (E — 1— 2k) ax) y* = 0. 
k=0 


We have thus obtained the recursion relation 


teig e LA2N-E 
ine (EEO Ve 


k =0,1,2,.... (8.76) 


Now recall that apart from a possible multiplicative polynomial factor, u(y) ~ ey /2 
or ew" /2, Therefore, H(y) either goes like a polynomial times eY ora polynomial. 
The function H(y) will be a polynomial only if € takes on specific values which 
truncate the infinite power series. 

Before we jump to this conclusion, I want to show that if the power series does 
not truncate, then H(y) does indeed grow like ev” for large y. We first write the 
Taylor series for ey: 


=) ta, 
k=0 


where 
1 


b= 2D! neven; 
k= 
0, kodd. 
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For large values of y, the later terms of the power series dominate the earlier terms 
and the behavior of the function is determined by the large k coefficients. Provided 
k is even, we have 


be+2 _ __(k/2)! 
be ((k+2)/2)! 
1 
~ (k/2) +1 
2 
~E. 


The ratio of the coefficients would have been the same even if we had multiplied 
ey by a polynomial. 

The behavior of the coefficients for large k of our solution to the Schrédinger 
equation is given by the recursion relation (8.76, and is 


Ar+2 _ (1 + 2k) -—€E 
am  (k+2)(kK+1) 
2k 
~ Re 
_ 2 
=e 


the same behavior as the Taylor series for e!”. For large y, then, the infinite power 
series for H(y) will grow as eY” times a polynomial. Since this does not satisfy the 
boundary conditions for u(y) at infinity, we must force the power series to truncate, 
which it does if we set € = €,, where 


Ey =14+2n, n=0,1,2,..., 


resulting in quantization of the energy. The dimensional result is 


with the ground state energy level given by Ey = fiw/2. The wavefunctions asso- 
ciated with each E, can be determined from the power series and are the so-called 
Hermite polynomials times the decaying exponential factor. The constant coefficient 
can be determined by requiring the wavefunctions to integrate to one. 

For illustration, the first two energy eigenfunctions, corresponding to the ground 
state and the first excited state, are given by 


mw\1/4 _ mw\1/4 /2mw  _, 442 
(x) _ (=) eT mux [mh (x) = ( ~) " xe mux [2h 


8.8.5 Particle in a three-dimensional box 


To warm up to the analytical solution of the hydrogen atom, we solve what may be 
the simplest three-dimensional problem: a particle of mass m able to move freely 
inside a cube. Here, with three spatial dimensions, the potential is given by 


0, O<x,y,z<L, 
co, otherwise. 


Vixy,2) = 
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We may simply impose the boundary conditions 


$(0,y,2) = PL, ¥,2) = p(x,0,2) = p(x, L,z) = p(x y,0) = p(x y,L) = 0. 


The Schrédinger equation for the particle inside the cube is given by 


rr (ep ep ey ole 
2m \ ox2 | oy? ' Qz2 Y. 


We separate this equation by writing 


p(x, y,Z) = X(x)¥(y)Z(z), 


and obtain 
2 


hi 
5 (X"YZ + XY"Z + XYZ") = EXYZ. 


Dividing by XYZ and isolating the x-dependence first, we obtain 


h x” i yl VAL 
2m X 2m ( Yo Z ) 
The left-hand side is independent of y and z and the right-hand side is independent 
of x so that both sides must be a constant, which we call E;. Next isolating the 
y-dependence, we obtain 


a y" i VAL 
— t+ E—E,. 
2m ( Y ) 2m ( Z ) 
The left-hand side is independent of x and z and the right-hand side is independent 


of x and y so that both sides must be a constant, which we call Ey. Finally, we define 
E, = E—E, — Ey. The resulting three differential equations are given by 


2mE, 
i 


eis Z=0. 
h 


2mE 
SH; ee el vie 


These are just three independent one-dimensional box equations so that the energy 
eigenvalue is given by 


(n2 + n2 + 02) nh? 
2mL? 


Enynynz 


and the associated wavefunction is given by 


3/2. i : 
~ (2) sin “ sin a sin", 0< x,y,z < L; 
NyNyn : 
vege 0, otherwise. 


8.8.6 The hydrogen atom 


Hydrogen-like atoms, made up of a single electron and a nucleus, are the atomic 
two-body problems. As is also true for the classical two-body problem, consisting 
of say a planet and the Sun, the atomic two-body problem can be reduced to a 
one-body problem by transforming to center-of-mass coordinates and defining a 
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Figure 8.5: The spherical coordinate system, with radial distance r, polar angle @ and 
azimuthal angle ¢. 


reduced mass pt. We will not go into these details here, but will just take as the 
relevant Schrodinger equation 


— 
a ra p+V(r)p = Ey, (8.77) 


where the potential energy V = V(r) is a function only of the distance r of the 
reduced mass to the center-of-mass. The explicit form of the potential energy from 
the electrostatic force between an electron of charge —e and a nucleus of charge 
+Ze is given by 

Ze 
~ Artegr’ 


V(r) = (8.78) 


With V = V(r), the Schrédinger equation (8.77) is separable in spherical coordi- 
nates. With reference to Fig. 8.5, the radial distance r, polar angle @ and azimuthal 
angle ¢ are related to the usual cartesian coordinates by 


x=rsinOcos?, y=rsin@sing, z—=rcosé; 
and by a change-of-coordinates calculation, the Laplacian can be shown to take the 


form 
ae 1 28 a 1 9 
) a 2 . 
V r2 or (r = T 7) sin @ 00 (sino) T ) ane r ag : (8.79) 


The volume differential dt in spherical coordinates is given by 


dt = r* sin Odrdodo. (8.80) 
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A complete solution to the hydrogen atom is somewhat involved, but neverthe- 
less is such and important and fundamental problem that I will pursue it here. Our 
final result will lead us to obtain the three well-known quantum numbers of the 
hydrogen atom, namely the principle quantum number n, the azimuthal quantum 
number J, and the magnetic quantum number m. 

With wp = (1,0, ¢), we first separate out the angular dependence of the Schrédinger 
equation by writing 


P(r, 8,p) = R(r)Y(6,¢). (8.81) 


Substitution of (8.81) into (8.77) and using the spherical coordinate form for the 
Laplacian (8.79) results in 


RY d (paR\ | RO. ay) | oR #Y 
Qu [12 dr dr) ' r2sin000 00) | 72sin2 9 ag? 
+ V(r)RY = ERY. 


To finish the separation step, we multiply by —2 yr? /it” RY and isolate the r-dependence 
on the left-hand side: 


1[d (dR 2Qur? 1 0. Oy 2.20 L. oy 
a 0 : 
R E (« dr )| gp YO) ="y | anand ad) ant oe 
The left-hand side is independent of 6 and ¢ and the right-hand side is independent 


of r, so that both sides equal a constant, which we will call A,. The R equation is 
then obtained after multiplication by R/ 7: 


1d (dR 2u Ay 
R=0. 82 
v2 dr (" 7) \is as) 72 (882) 
The Y equation is obtained after multiplication by Y: 
1 a oY 1 #y 
in @ + A,Y =0. : 
sind 00 (sin a0 ) neue ee) 


To further separate the Y equation, we write 
Y(8,4) = (8) ®(4). (8.84) 


Substitution of (8.84) into (8.83) results in 


® d dO © do 
ind + A,;O = 0. 
sind do (: 7) ae lala 
To finish this separation, we multiply by sin? @/@® and isolate the 6-dependence 
on the left-hand side to obtain 


sndd /(.. dO\ , .5, 1d 
© do (sino ) + Ay sin 9= 8 age 


The left-hand side is independent of ¢ and the right-hand side is independent of 
8 so that both sides equal a constant, which we will call Ay. The © equation is 
obtained after multiplication by ©/ sin 6: 


1 d/. dO\ | Ae 5 
aa (sino? ) (a {3,)e-0 (8.85) 
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The ® equation is obtained after multiplication by ®: 


a 
ag? +A2® = 0. (8.86) 
The three eigenvalue ode equations for R(r), O(@) and ®(#) are thus given by 
(8.82), (8.85) and (8.86), with eigenvalues E, A; and Ap. 
Boundary conditions on the wavefunction determine the allowed values for the 
eigenvalues. We first solve (8.86). The relevant boundary condition on ® = ®(¢) is 
its single valuedness, and since the azimuthal angle is a periodic variable, we have 


&(p +27) = (p). (8.87) 


Periodic solutions for ®(@) are possible only if Az > 0. We therefore obtain using 
the complex form the general solutions of (8.86): 
AeiV429 + Be-iVA2h, Ay > 0, 
O(p) = a 
C+D¢, Ay = 0. 
The periodic boundary conditions given by (8.87) requires that \/A2 is an integer 


and D = 0. We therefore define Az = m2, where m is any integer, and take as our 


eigenfunction 


1 3 
Di, (p) = am” 


where we have nomalized ®,, so that 


27 
| |®n(p) Pd = 1. 


The quantum number m is commonly called the magnetic quantum number be- 
cause when the atom is placed in an external magnetic field, its energy levels be- 
come dependent on m. 

To solve the © equation (8.85), we let 


w=cosé, P(w) = (6). 


Then 
sin? @ =1—w* 
and 
d® dPdw . ,dP 
i eae a 
allowing us the replacement 
: d 
7 a sin d=. 


With these substitutions and A7 = m7”, (8.85) becomes 
dP 


d 2 m2 
ce a w Ce (a 2) P=0. (8.88) 


To solve (8.88), we first consider the case m = 0. Expanding the derivative, (8.88) 


then becomes 4 
d“P dP 
(1—@) aa 20, tA,P=0. (8.89) 
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Since this is an ode with nonconstant coefficients, we try a power series ansatz of 
the usual form 


P(w)=)° awk. (8.90) 
k=0 
Substitution of (8.90) into (8.89) yields 
y k(k = 1)agwk* — SO k(k — 1) age — YP 2ka,w* + D> Ayazw* = 0. 
k=2 k=0 k=0 k=0 


Shifting the index in the first expression and combining terms results in 


foe) 


Y {tk + 2)(k + 1)aps2 — [(k(k +1) Ajay fo = 0. 


k=0 


Finally, setting the coefficients of the power series equal to zero results in the recur- 
sion relation 
P _ k(k+1)-”4 7 
ota (so) ee aa 


Even and odd coefficients decouple, and the even solution is of the form 
CO 
Peven(w) = ye Arnw", 
n=0 
and the odd solution is of the form 
= 1 
Poga(w) = > dong wrt? 


n=0 


An application of Gauss’s test for series convergence can show that both the even 
and the odd solutions diverge when |w| = 1 unless the series terminates. Since the 
wavefunction must be everywhere finite, the series must terminate and we obtain 
the discrete eigenvalues 


A, =1(1+1), for! =0,1,2,.... (8.91) 


The quantum number / is commonly called the azimuthal quantum number de- 
spite having arisen from the polar angle equation. The resulting eigenfunctions 
P;(w) are called the Legendre polynomials. These polynomials are usually normal- 
ized such that P;(1) = 1, and the first four Legendre polynomials are given by 


Po(w) = 1, Pi(w) =w, 
P,(w) = 5 (30? —1), P3(w) = 5 (50? — 30). 


With A; = 1(1 +1), we now reconsider (8.88). Expanding the derivative gives 


ar dP m 
2 
=0. 92 
(l-w a 2we (ww 1) 4) P=0 (8.92) 
Equation 8.92) is called the associated Legendre equation. The associated Legendre 
equation with m = 0, 

d?P dP 


(1 wT 207 +11 +1)P =0. (8.93) 
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is called the Legendre equation. We now know that the Legendre equation has 
eigenfunctions given by the Legendre polynomials, P;(w). Amazingly, the eigen- 
functions of the associated Legendre equation can be obtained directly from the 
Legendre polynomials. For ease of notation, we will assume that m > 0. To include 
the cases m < 0, we need only replace m everywhere by |m|. 

To see how to obtain the eigenfunctions of the associated Legendre equation, we 
will first show how to derive the associated Legendre equation from the Legendre 
equation. We will need to differentiate the Legendre equation (8.93) m times and to 
do this we will make use of Leibnitz’s formula for the mth derivative of a product: 


Sede = ¥ (") Ee 


m= \i dxi dx™—i" 


where the binomial coefficients are given by 


(“) _ m! 
j fida—pr 
We first compute 


qm »@P] et (m\ [ di >») [ ai ap 
dw™ a . ed » (") Ee vam) dw™—i dw? | ° 


Only the terms j = 0,1 and 2 contribute, and using 


(=e Gm GS 


and the more compact notation 


a"P _ pin) 
dw" p(w), 
we find 
_ l(a = w?) P| = (1— w*) PO") — amwP" — m(m—1)P™. — (8.94) 


We next compute 


a™ —. dP) & (m\ [ di d”-i dP 
iow Pe] = 25 (7) La Lior 


Here, only the terms j = 0 and 1 contribute, and we find 
ii (1) (m+1) (m) 
= m m 
= [2wP 2wP(™+1) 4 9mp(™), (8.95) 
Differentiating the Legendre equation (8.93) m times, using (8.94) and (8.95), and 
defining 
p(w) = P\” (w) 
results in 
dp 
dw 


2(1-+ m)w 2? + (1-41) —m(m-+ 1] p =0. (8.96) 


(1 w’) a 
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Finally, we define 


q(w) = (1— w?)"/? p(w). (8.97) 
Then using 
ap 4 ay—my2 ( 49 mwq 
ane) aa Gaae (8.98) 
and 


dp 2\—-m/2 fq , 2mw dq , [m(m+2)w? — m 
Fr al \ aad '1—w dw | | (1—w)2 | ala} wee) 


we substitute (8.98) and (8.99) into (8.96) and cancel the common factor of (1 — 
w?)—™/2 to obtain 


a m? 
hie = 
(l1-w ) Fo 2we fu + 1) | q=0, 


which is just the associated Legendre equation (8.92) for q = q(w). We may call 
the eigenfunctions of the associated Legendre equation P,,,(w), and with q(w) = 
Pin (w), we have determined the following relationship between the eigenfunctions 
of the associated Legendre equation and the Legendre polynomials: 


P =(1—-w* miso 
1m (@) = ( w ) d 
w 


Pi(w), (8.100) 


where we have now replaced m by its absolute value to include the possibility of 
negative integers. Since P;(w) is a polynomial of order I, the expression given by 
(8.100) is nonzero only when |m| < 1. Sometimes the magnetic quantum number m 
is written as m, to signify that its range of allowed values depends on the value of 
L. 

At last, we need to solve the eigenvalue ode for R = R(r) given by (8.82), with 
V(r) given by (8.78) and A, given by (8.91). The radial equation is now 


Ld f oA Qu Ze* I(1+1) 
E+ R=0. 101 
r? dr (r dr ) { i | Arter 72 : ert) 
Note that each term in this equation has units of one over length squared times R 
and that E < 0 for a bound state solution. 
It is customary to nondimensionalize the length scale so that 2vE/h? = —1/4 in 


dimensionless units. Furthermore, multiplication of R(r) by r can also simplify the 
derivative term. To these ends, we change variables to 


_ VS 8ulE| 


p= Er, ule) = rR(0), 


and obtain after multiplication of the entire equation by p the simplified equation 


ee, [(e 1. Ue). 
we {3 oar humo, (8.102) 


where a now plays the role of a dimensionless eigenvalue, and is given by 


Ze* il 
«= Fon DE): (8.103) 
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As we saw for the problem of the simple harmonic oscillator, it may be helpful 
to consider the behavior of u = u(p) for large p. For large p, (8.102) simplifies to 


@u 1 
dp? a aa 0, 


with two independent solutions 
ee /2, 
u(p) = oe 


Since the relevant boundary condition here must be limp-co u(p) = 0, only the 
second decaying exponential solution can be allowed. 

We can also consider the behavior of (8.102) for small p. Multiplying by 7, 
and neglecting terms proportional to p and 0" that are not balanced by derivatives, 
results in the Cauchy-Euler equation 


au 
2 _ 
p ae 1(i+1)u=0, 


which can be solved by the ansatz u = p*. After canceling 0°, we obtain 


s(s—1)—-1(1+1) =0, 


1+1, 
— 
—l. 


If R = R(r) is finite at r = 0, the relevant boundary condition here must be 
limp_,9 u(e) = 0 so that only the first solution u(o) ~ o'*? can be allowed. 

Combining these asymptotic results for large and small p, we now try to substi- 
tute 


which has the two solutions 


u(p) = p!*1e-P/?F() (8.104) 
into (8.102). After some algebra, the resulting differential equation for F = F(p) is 
found to be ; 

ee | (74 1) dF | a—~U+1)p_o9 
do p dp p 


We are now in a position to try a power-series ansatz of the form 


F(p) _ z axp* 
k=0 


to obtain 
k(k — 1)a,o*-? + > 2(1+ 1)ka,o*-? _ > kayo + y ja —(1+1)] a,pk—! =0. 
k=2 k=1 k=1 k=0 


Shifting indices, bringing the lower summations down to zero by including zero 
terms, and finally combining terms, we obtain 


yf [(k + 1)(k+2(1+1)) Jays —(14+!4+k- wa boh =0. 
k=0 
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Setting the coefficients of this power series equal to zero gives us the recursion 


relation 
1+l+k—-«a 


“e+ ep a)(R+200+1) - 


For large k, we have ax+1/a, —> 1/k, which has the same behavior as the power 
series for e?, resulting in a solution for u = u(o) that behaves as u(o) = e?/? for 
large p. To exclude this solution, we must require the power series to terminate, 
and we obtain the discrete eigenvalues 


k- 


a=1+4+l+4+n', n'’=0,1,2,.... 


The function F = F(p) is then a polynomial of degree n’ and is known as an asso- 
ciated Laguerre polynomial. 

The energy levels of the hydrogen-like atoms are determined from the allowed 
« eigenvalues. Using (8.103), and defining 


n=1+l+4mn’, 
for the nonnegative integer values of n' and I, we have 
uZet 


En = sy 2 =1,2,3,.... 
2(A7teg)*h*n2 


If we consider a specific energy level E,,, then the allowed values of the quantum 
number ! are nonnegative and satisfy 1 = n — n’ — 1. For fixed n then, the quantum 
number | can range from 0 (when n’ = n — 1) ton — 1 (when n’ = 0). 

To summarize, there are three integer quantum numbers n, 1, and m, with 


n= 12,3524; 
1=0,1,...,n—1, 
m=-l,...,l, 
and for each choice of quantum numbers (n,1,m) there is a corresponding energy 
eigenvalue E,, which depends only on n, and a corresponding energy eigenfunction 
Y = Wnim(t,9,), which depends on all three quantum numbers. 
For illustration, we exhibit the wavefunctions of the ground state and the first 


excited states of hydrogen-like atoms. Making use of the volume differential (8.80), 
normalization is such that 


co pm p27 
7 | i lWrtm (1, 9, )|71? sin 0 dr dé dp = 1. 
0 0 JO 


Using the definition of the Bohr radius as 


Amteoh* 
ag = 5] 
pe 


, 


the ground state wavefunction is given by 


a Pag 
_ —Zr/ag 
Ses e : 
ve Vr (=) 


and the three degenerate first excited states are given by 


i yee Zr\ _7 
Se 7 eran er 712A 
caaae (x) ( 
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